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(57)Abstract: 

PROBLEM TO BE SOLVED: To control a coupling direction of a 
ferroelectric insulation film, restrict mutual operation with an 
electric field and restrict a deterioration in an induced 
electricity due to the electric field, by a method wherein 
simulation resultants, based on molecular dynamics and fluid 
models, are utilized for other simulations. 
SOLUTION: In a process fluctuation calculation unit I, macro- 
property not taking into consideration fundamentally an atomic 
level is calculated, and is stored in a file 2. In file contents, a 
coupling angle of an interatomic position vector calculated in a 
molecular dynamics simulation II is inputted to a film 
characteristic calculation unit III, and a fundamental film 
characteristic is calculated therein and returned to a file 3. The 
fundamental film characteristic calculated in the film 
characteristic calculation unit III is also returned to the process 
fluctuation calculation unit I as a part of parameters. The 
molecular dynamics simulation II is executed in combination 
with a molecular dynamics method and a density function. As a 
result, it is possible to progress a semiconductor manufacturing 
process according to a desirable process, or with amendments without a test piece. 
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CLAIMS 

[Claim(s)] 

[Claim 1]In a simulation method of physical development in a manufacturing installation ot a 
manufacturing process of electronic parts, A manufacturing process factor about the characteristic of 
an electronic industry material formed in said manufacturing installation is made into an input value, Are 
prediction data of the characteristic of an electronic industry material which performed a three- 
dimensional simulation based on this input value, and was formed in a manufacturing installation the 
method of obtaining, and said simulation, It consists of a process performed based on molecular 
dynamics, and a process performed based on a fluid model, At least one side of three-dimensional 
physical quantity outputted as a simulation result performed based on said molecular dynamics and 
three-dimensional physical quantity outputted as a simulation result performed based on said fluid 
model, A designing method of a manufacturing process of electronic parts which being transmitted to 
another side and using for a simulation of a there. 

[Claim 2]A designing method of the electronic-parts manufacturing process according to claim 1 
characterized by transmitting physical quantity, such as a viscoelasticity coefficient, a Poisson's ratio, a 
diffusing constant, and Cij, from said three-dimensional molecular dynamics simulation to a process 
simulation of said three-dimensional fluid model. 

[Claim 3] A designing method of the electronic-parts manufacturing process according to claim 1 
characterized by transmitting physical quantity, such as a stress tensor, temperature, impurity 
concentration, and point defect concentration, from a process simulation of said three-dimensional fluid 
model to said three-dimensional molecular dynamics simulation. 

[Claim 4]A designing method of the electronic-parts manufacturing process according to claim 1, 
wherein said electronic parts are semiconductor devices. 

[Claim 5]A designing method of the electronic-parts manufacturing process according to claim 1, 
wherein it has a stage of inputting fluctuation of said preset value and said molecular dynamics 
simulation performs variations calculation based on said input value and its fluctuation. 
[Claim 6]A designing method of the electronic-parts manufacturing process according to claim 1, 
wherein said molecular dynamics simulation is performed in combination of a molecular dynamics 
method and a density functional method. 

[Claim 7]A designing method of the electronic-parts manufacturing process according to claim 6, 
wherein said molecular dynamics method is based on the Takashina round partial differential method. 
[Claim 8]A designing method of the electronic-parts manufacturing process according to claim 7, 
wherein said Takashina round partial differential method is dealing with a three-body problem. 
[Claim 9]A designing method of the electronic-parts manufacturing process according to claim 1, 
wherein said molecular dynamics simulation is performed using experiential potential. 
[Claim 10]A designing method of the electronic-parts manufacturing process according to claim 9, 
wherein said experiential potential is ATTA potential. 

[Claim 11]A designing method of the electronic-parts manufacturing process according to claim 9, 
wherein said experiential potential is BKS potential. 

[Claim 12]A process of measuring the characteristic of an electronic industry material formed in a 



manufacturing installation in a manufacturing process of electronic parts, and obtaining real 
observational data, A process of obtaining prediction data of the characteristic of an electronic industry 
material which made an input value a manufacturing process factor about the characteristic of an 
electronic industry material formed in said manufacturing installation, performed a three-dimensional 
simulation based on this input value, and was formed in said manufacturing installation, By a process of 
carrying out comparison assay of said prediction data and the real observational data, and said 
comparison assay. When a significant difference is accepted between a preset value of said 
manufacturing process factor, and said manufacturing process factor guessed from said real 
observational data, Have a process of carrying out correction processing of said manufacturing process 
factor, and said simulation, It consists of a process performed based on molecular dynamics, and a 
process performed based on a fluid model, At least one side of three-dimensional physical quantity 
outputted as a simulation result performed based on said molecular dynamics and three-dimensional 
physical quantity outputted as a simulation result performed based on said fluid model, A manufacturing 
method of electronic parts which being transmitted to another side and using for a simulation of a there. 

[Claim 13]A manufacturing method of the electronic parts according to claim 12 a significant difference 
between said manufacturing process factors guessed from said real observational data having crossed 
tolerance level, and processing as inferior goods when correction processing is impossible. 
[Claim 14]A manufacturing method of the electronic parts according to claim 12 currently carrying out 
clustering of said manufacturing installation. 

[Claim 15]A manufacturing method of the electronic parts according to claim 14, wherein said 
manufacturing installation which carried out clustering is single wafer processing. 
[Claim 16]A manufacturing method of the electronic parts according to claim 12 characterized by 
transmitting physical quantity, such as a viscoelasticity coefficient, a Poisson's ratio, a diffusing 
constant, and Cij, from said three-dimensional molecular dynamics simulation to a process simulation of 
said three-dimensional fluid model. 

[Claim 17]A manufacturing method of the electronic parts according to claim 12 characterized by 
transmitting physical quantity, such as a stress tensor, temperature, impurity concentration, and point 
defect concentration, from a process simulation of said three-dimensional fluid model to said three- 
dimensional molecular dynamics simulation. 

[Claim 18]A manufacturing method of the electronic parts according to claim 12, wherein said 
characteristics are a film formed on a semiconductor substrate, a semiconductor substrate surface, and 
the characteristic about at least one inside a semiconductor substrate. 

[Claim 19]A manufacturing method of the electronic parts according to claim 12, wherein comparison 
assay with said prediction data and real observational data is performed in real time one by one. 
[Claim 20]A manufacturing method of the electronic parts according to claim 12, wherein a process of 
carrying out correction processing of said manufacturing process factor is performed in real time one by 
one. 

[Claim 21]A manufacturing method of the electronic parts according to claim 12, wherein said molecular 
dynamics simulation is performed in combination of a molecular dynamics method and a density 
functional method. 

[Claim 22]A manufacturing method of the electronic parts according to claim 21, wherein said molecular 
dynamics method is based on the Takashina round partial differential method. 

[Claim 23]A manufacturing method of the electronic parts according to claim 22, wherein said Takashina 
round partial differential method is dealing with a three-body problem. 

[Claim 24]A manufacturing method of the electronic parts according to claim 12, wherein said molecular 
dynamics simulation is performed using experiential potential. 

[Claim 25]A manufacturing method of the electronic parts according to claim 24 characterized by a 
thing in which said experiential potential is ATTA potential, and to do. 

[Claim 26]A manufacturing method of the electronic parts according to claim 24 characterized by a 
thing in which said experiential potential is BKS potential, and to do. 

[Claim 27]A manufacturing method of the electronic parts according to claim 12 characterized by 
performing a process simulation of a fluid model using at least one prediction data of two or more of 
said characteristics obtained in said molecular dynamics simulation. 



[Claim 28]A manufacturing method of the electronic parts according to claim 27, wherein both said 
molecular dynamics simulation and a process simulation of a fluid model are three-dimensional 
simulations. 

[Claim 29]A manufacturing method of the electronic parts according to claim 27 characterized by 
transmitting physical quantity, such as a viscoelasticity coefficient, a Poisson's ratio, a diffusing 
constant, and Cij, from said three-dimensional molecular dynamics simulation to a process simulation of 
said three-dimensional fluid model. 

[Claim 30]Have the following, and said computer system performs said prediction data and real 
observational data, and comparison assay by said comparison assay. A manufacturing system of 
electronic parts characterized by carrying out correction processing of the manufacturing process 
factor via a control means of said manufacturing installation when a significant difference is accepted 
between a preset value of said manufacturing process factor, and said manufacturing process factor 
guessed from said real observational data. 
A manufacturing installation of electronic parts. 

A measuring device which measures the characteristic of an electronic industry material of said 
electronic parts formed in said manufacturing installation. 

A computer system which makes an input value a preset value of a manufacturing process factor of 
said electronic parts, performs a simulation based on this input value, and obtains prediction data of 
said characteristic. 

A means to connect said computer system with said measuring device and said manufacturing 
installation. 

[Claim 31]A manufacturing system of the electronic parts according to claim 30 currently carrying out 
clustering of said manufacturing installation. 

[Claim 32]A manufacturing system of the electronic parts according to claim 31, wherein said 
manufacturing installation which carried out clustering is single wafer processing. 

[Claim 33]A manufacturing system of the electronic parts according to claim 30, wherein said simulation 
is performed based on a molecular dynamics simulation by combination of a molecular dynamics method 
and a density functional method. 

[Claim 34]A manufacturing system of the electronic parts according to claim 30, wherein said 
characteristics are a film formed on a semiconductor substrate, a semiconductor substrate surface, and 
the characteristic about at least one inside a semiconductor substrate. 

[Claim 35]A manufacturing system of the electronic parts according to claim 30, wherein comparison 
assay with said prediction data and real observational data is performed in real time one by one. 
[Claim 36]A manufacturing system of the electronic parts according to claim 30, wherein a process of 
carrying out correction processing of said manufacturing process factor is performed in real time one by 
one. 

[Claim 37]A manufacturing system of the electronic parts according to claim 30, wherein fluctuation of 
said preset value is inputted as said input value and said molecular dynamics simulation performs 
variations calculation based on said input value and its fluctuation. 

[Claim 38]A manufacturing system of the electronic parts according to claim 30, wherein said molecular 
dynamics method is based on the Takashina round partial differential method. 

[Claim 39]A manufacturing system of the electronic parts according to claim 38, wherein said Takashina 
round partial differential method is dealing with a three-body problem. 

[Claim 40]A manufacturing system of the electronic parts according to claim 30, wherein said molecular 
dynamics simulation uses together a method which used experiential potential. 
[Claim 41 ]A manufacturing system of the electronic parts according to claim 40 characterized by a 
thing in which said experiential potential is ATTA potential, and to do. 

[Claim 42]A manufacturing system of the electronic parts according to claim 40 characterized by a 
thing in which said experiential potential is BKS potential, and to do. 

[Claim 43]A manufacturing system of the electronic parts according to claim 30 characterized by 
performing a process simulation of a fluid model using at least one prediction data of said characteristic 
obtained in said molecular dynamics simulation. 

[Claim 44]A manufacturing system of the electronic parts according to claim 43, wherein both said 



molecular dynamics simulation and a process simulation of a fluid model are three-dimensional 
simulations. 

[Claim 45]A manufacturing system of the electronic parts according to claim 43 characterized by 
transmitting physical quantity, such as a viscoelasticity coefficient, a Poisson's ratio, a diffusing 
constant, and Cij, from said three-dimensional molecular dynamics simulation to a process simulation of 
said three-dimensional fluid model. 

[Claim 46]A procedure of inputting a preset value of a manufacturing process factor about the 
characteristic of an electronic industry material which is a program which performs a simulation of 
physical development in a manufacturing installation of a manufacturing process of electronic parts, and 
was formed in said manufacturing installation, A recording medium which recorded a program which 
performs a procedure of performing a molecular dynamics simulation by combination of a molecular 
dynamics method and a density functional method based on this input value, and obtaining prediction 
data of said characteristic and in which computer reading is possible. 

[Claim 47]The recording medium according to claim 46, wherein said molecular dynamics simulation is 
based on combination of a molecular dynamics method and a density functional method. 
[Claim 48]The recording medium according to claim 47, wherein said molecular dynamics method is 
based on the Takashina round partial differential method. 

[Claim 49]The recording medium according to claim 48, wherein said Takashina round partial differential 
method is dealing with a three-body problem. 

[Claim 50]The recording medium according to claim 47, wherein said electronic parts are semiconductor 
devices. 

[Claim 51]The recording medium according to claim 47, wherein it has a procedure of inputting 
fluctuation of said preset value and said molecular dynamics simulation performs variations calculation 
based on said input value and its fluctuation. 

[Claim 52]The recording medium according to claim 47, wherein said molecular dynamics simulation uses 
together a method which used experiential potential. 

[Claim 53]The recording medium according to claim 52 characterized by a thing in which said 
experiential potential is ATTA potential, and to do. 

[Claim 54]The recording medium according to claim 52 characterized by a thing in which said 
experiential potential is BKS potential, and to do. 

[Claim 55]The recording medium according to claim 47 characterized by performing a process simulation 
of a fluid model using prediction data of said characteristic obtained in said molecular dynamics 
simulation. 

[Claim 56]The recording medium according to claim 55, wherein both said molecular dynamics simulation 
and a process simulation of a fluid model are three-dimensional simulations. 

[Claim 57]The recording medium according to claim 55 characterized by transmitting physical quantity, 
such as a viscoelasticity coefficient, a Poisson's ratio, a diffusing constant, and Cij, from said three- 
dimensional molecular dynamics simulation to a process simulation of said three-dimensional fluid 
model. 
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DETAILED DESCRIPTION 



[Detailed Description of the Invention] 
[0001] 

[Field of the Invention]This invention relates to the manufacturing method of electronic parts, a 
manufacturing system, a designing method, and a recording medium, A limited quantity in single wafer 
processing or the semiconductor device manufacturing installation which carried out clustering, or the 
"spot" measured value of qualitative contents is especially extended to the maximum, A semiconductor 
device manufacturing process is related with the manufacturing method of the semiconductor device 
which makes it possible to go on correcting, the manufacturing system, designing method, and recording 
medium as a desired process, without using a test piece. 
[0002] 

[Description of the Prior Art]Conventionally, the manufacturing installation of the semiconductor device 
has continued up to now through the following changes. If this is looked at in each generation's flow, the 
generation of dRAM and the chip cost for every generation of the will change, as shown in Drawing 128, 
while capacity increases from 64K to 256M. The vertical axis of Drawing 128 is expressing the chip cost 
in 1M as 1. In Drawing 128, the plotted point shows chip cost and the horizontal axis shows the 
generation of dRAM. 

[0003]As shown in Drawing 128, also as for the routing counter, 64K to 256K is comparatively stable, 
and its chip cost is also stable. Whenever it passes through a generation from 4M, the routing counter is 
increasing and chip cost is also going up in connection with it. In order to reduce the increase in cost by 
increase of this routing counter, as shown also in Drawing 128, in 64M. (1) Means, such as 
standardization [ of useless exclusion and (2) work ], reduction [ of (3) chip sizes ], (4) minuteness- 
making promotion, (5) productivity drive, and (6) FA-ized promotion and formation of (7) large calibers, 
were provided, and cost reduction is measured. 

[0004]If set to 256M, it will put on "Bai" Rule by attaining application of an exotic material, adoption of 
a high dielectric, application of new structure, solidification, etc. However, if it becomes 1G class, a 
situation will turn around and will become [ whether "Bai" Rule is followed and ] unknown. In the 
conventional method of going in pursuit of highly-efficient-izing of manufacturing installation each to 
256M old class to a limit, it becomes difficult to already calculate the solution as an enterprise. It will 
cause going in pursuit of highly-efficient-izing of manufacturing installation each to a limit, i.e., a high 
cost. 

[0005]As shown in Drawing 129, for every process, from the process A, the device made highly efficient 
shall be used in order of B->C->D->E->F, and shall go, for example. In this case, it is clear that 
increase of an equipment field also causes a high cost. 

[0006]His not forgetting in more than one more 1G class is that the close control by physical necessity 
becomes indispensable in each process. For example, they are interception of contamination, exclusion 
of a natural oxidation film, pursuit of flattening in an atom level, etc. In the present, in consideration of 
these things, although it is a research stage, the manufacturing installation which carried out clustering 
as shown in Drawing 130 is proposed as a new concept. This manufacturing installation (cluster tool) 
that carried out clustering is mainly single wafer processing, although there is also a batch type. 



As shown also in Drawing 130, there is a conveyance system in the center, like B->C->D from A, it 
goes on one by one and each process goes, respectively, as shown in an arrow. 
In many cases, it is maintained at a high vacuum, and covering of a natural oxidation film can also 
decrease, and it has come to be able to perform interception of contamination, and surface control with 
an atom level within each process and each process at this time. 

[0007]On the other hand, vigor and a ferroelectric are going to thin film-ization about the electric field 
effect MIS element formed on the semiconductor in the technical development of improvement in the 
speed of these days, and high integration, natural, when a ferroelectrics film is thin-filrrHzed 
threshold voltage becomes shallow, and main movement speed boils especially AC characteristics 
markedly by this early, and is made to improve only at this rate If considering EEPROM etc. minuteness 
making is carried out, it will become a very severe operating condition. In such a case, in the 
ferroelectrics film of the conventional usual process, sufficient reliability is not acquired any longer. 
[0008]However, even if minuteness making of the MIS element is carried out, the remedy of the 
characteristic of the ferroelectrics film itself is seldom considered, but spends for a high electric field to 
a gate ferroelectrics film especially in the state where power supply voltage does not fall so much, at 
the time of main movement. The electron hole generated in impact ionization etc. from channel regions 
is poured in into a ferroelectrics film by boundary conditions, such as the polarity of a gate electrode, 
and drain voltage, respectively. And into a ferroelectrics film, the trap of these careers is carried out 
and it not only reduces long term reliability, but they cause a resisting pressure fall etc. 
[0009]It sees with an atom level to a slight degree. In a heat ferroelectrics film, if a high electric field is 
impressed, for example to a silicon ferroelectrics film, Si-0 combination which constitutes a network will 
interact with the high electric field impressed from the outside. And combination is cut, the trapping 
center where an electron and an electron hole are captured is formed, it is captured in this trapping 
center, the field strength distribution of a thickness direction becomes higher than an average electric 
field locally, and the electron and electron hole through which it passes continuously result in a 
dielectric breakdown soon. 

[0010]lt is going to improve these situations, and these days, although the idea which single-crystal- 
izes a ferroelectrics film was a society level, it began to be advocated. For example, it is reported to 
J.Appl.Phys., Aol.69 (12), and p8313 (1991) that cerium oxide Ce0 2 carries out single crystal growth on 
a silicon (111) side. Or it is reported to Japan. J.Appl.Phys.Suppl. and Aol.21-1 and p187 (1982) that 
calcium fluoride CaF2 carries out single crystal growth to silicon single crystal Kami. These are parts. 
It is not only still in the state of an idea, but the question remains in the calculation technique which 
becomes a basis of the indicator plentifully. 

Then, even if it raises a ferroelectrics film to an example, for example, it is in the situation where the 
figure of the single crystal of the ferroelectrics film using it itself is not recognized enough one. I will 
look at what was attached and reported to the structure of this single crystal ferroelectrics film. To 
one. A. Miyamot, K.Takeichi, T.Hattori, M.Kubo, and T.Inui, "Mechanism of Layer-Layer Homoepitaxial 
Growth of SrTiO3(100)as Investigated by Molecular Dynamics and Computer Graphics". Jpn. J. Appl. 
Phys. Vol 31 (1992), There is 4463-4464. Here, stable structure is calculated when pasting it up with a 
ground Si substrate by using perovskite as a single crystal ferroelectrics film. However, the angle 
between Ti-Si-Sr and the initial arrangement of distance are accidentally used for Si-O-Ti **. These 
have simplified the structure of a single crystal ferroelectrics film. 
The indicator is not necessarily right. 

[001 1]If it comes how to design a single crystal ferroelectrics film so that these things may show, the 
actual condition is almost right being unable to argue. The actual condition is still using the refinement 
on a formation process from beginning to end about improvement in the characteristic about another 
side and the ferroelectrics film itself. For example, the purest possible field is prepared beforehand or it 
is made the argument on the weld slag temperature for film formation, or sputtering ambient from 
beginning to end. It must be said that these have the still very low recognition including the structure 
about a single crystal ferroelectrics film or a ferroelectrics film. Since neither the Reason for the 
spontaneous polarization manifestation of a ferroelectric nor relation with the crystal structure is 
known, it passes only to phenomenological pursuit. 



[0012] 

[Problem(s) to be Solved by the Invention]However, in the above new cluster tools, the biggest difficulty 
faced now is the situation observation technique. That is, in the case of the device which performs 
single distance, the wafer first called a test piece was introduced into each process. Using this test 
piece, thickness, a refractive index, etc. were measured and it was being tested whether the process 
would advance as a predetermined condition. 

[0013]There is a problem on quantitative and qualitative understanding in monitoring of the conditions in 
each cluster tool by use of this test piece. For example, in the stage film formation or the growing 
process of an oxide film, usually uses the device of XPS or FT-IR. Though this is measured by in-situ, 
under the present circumstances, the place which the data means cannot fully be examined. It is the 
judgment grade whether the predetermined substance is generating at most judging from the peak wave 
number position of the spectrum. There are many places for which it depends on a researcher's 
judgment in such a situation. Quantitive correspondence cannot be performed and a still more prompt 
action cannot be performed, either. When there is a strange new substance for which it cannot depend 
on a researcher's experience, there is no specification of judgment. 

[0014]Use of the molecular dynamics simulator is also proposed by one side, however — comparing a 
molecular dynamics simulator with a process variation calculation program generally — (1) — 
calculation is dramatically slow and the execution area of the substance made into (2) objects is said to 
be very small. Then, most of a molecular dynamics simulator and the thing referred to as connecting a 
process variation calculation program directly was being unable to think. That is, it seemed that the 
molecular dynamics portion can pull the whole from speed with a late molecular dynamics portion after 
all, and is not fit for practical use even if it makes it unionHze since calculation is very slow. In a 
molecular dynamics portion, a calculation area is very small and is about at most tens of A, and on the 
other hand, since several micron field is calculated, it is difficult in a process variation calculation 
program, to connect these directly also from a point of a gap of a calculation area. 

[0015]For example, this invention persons were kept from forming level in an interface while this energy 
determined the position and made the highly reliable MIS element by this under impression conditions as 
installation requirements of a single crystal ferroelectrics film paying attention to the free energy of the 
whole system including an interface so that it might become the minimum. In searching for the energy of 
the whole system especially, although it needs to hold strictly and correctly reproducing the structure 
of a single crystal ferroelectrics film certainly, and how each of these constituted atoms exercise, This 
invention persons completed ab-initio molecular dynamics, further, in molecular dynamics, established 
the original potential estimation method, and enabled these calculations for the first time. And it is said 
that from a meaning to the degree of the quality of a single crystal ferroelectrics film which guarantees 
improving points, such as an electrical property, enough will be clarified using this. 

[0016]For example, combination which constitutes the basic structure of a tunnel insulating film if a high 
electric field is impressed to a tunnel insulating film, (for example, in the case of a silicon ferroelectrics 
film, it interacts with Si-0 joint), and combination is cut, the electron and electron hole which pass a 
tunnel insulating film following combination from which it was cut are captured, an electric field becomes 
larger than an average electric field locally according to the electric charge distribution in a film, and 
degradation advances further. In order to suppress advance of degradation, according to this invention 
person, it is suppressing cutting of combination (for example, Si-0 combination) of the insulator layer by 
this electric field. 

[001 7]It depends for the interaction with combination of an electric field and an insulator layer on 
direction of an electric field and the direction to combination of an insulator layer. Therefore, in order to 
weaken an interaction, the direction to combination of an insulator layer is making the thing of strong 
direction an interaction decrease more. So that amorphously, when completely random, a part of 
combination of an insulator layer serves as direction with a certainly strong interaction with an electric 
field. 

[0018]The purpose of this invention controls the direction to combination of a ferroelectric insulator 
layer, suppresses an interaction with an electric field, and aims at aiming at control of dielectric 
degradation by an electric field. 

[0019]This is related with high reliability dielectric film creation art, and this invention person etc. 
produce it based on the calculation result by the new abinitio molecular dynamics system, and original 



experimental data. The single crystal ferroelectrics film can get worse on the contrary than an 
amorphous dielectric film, if neither a substrate nor physical relationship with a substrate electrode 
besides selection of the material itself is optimized when this is applied to a ferroelectrics film etc. It is 
said that the material at this time and the criteria of crystal orientation will be shown. As a new 
indicator by this invention person, if a word dares to describe, it will be "performing crystal arrangement 
so that the free energy of a system may become the minimum under a main impressed electric field" 
etc. "Tolerance level etc. of the crystal defect concentration which may exist in a single crystal 
f erroe | ec trics film" is actually described as upper problem. "The tolerance level etc. where the free- 
energy distribution at the time of seeing locally is uneven in order that the free energy of the system at 
the time of main impression may deter the trigger of degradation other than the arrangement which 
becomes the minimum" is described. 

[0020]having very complicated structures, such as perovskite symmetrical with 4 times, even if it takes 
up even a single crystal ferroelectrics film, and the combination of a plane of composition with a 
substrate — large — the characteristic — things, although things etc. are very complicated, according 
to this invention person etc., "it gaining" under main movement voltage, even if "it loses" 
crystallographically — it comes out. 

[0021]This invention solves the problem of the monitoring facility in clustering, and A semiconductor 
device manufacturing installation, A limited quantity especially in single wafer processing or a clustering 
manufacturing installation or the "spot" measured value of qualitative contents is extended to the 
maximum, It aims at providing the manufacturing method of the semiconductor device which makes it 
possible to go on while a request corrects a semiconductor device manufacturing process without a 
test piece as a process, a manufacturing installation, a simulation method, and a simulator. 
[0022] 

[Means for Solving the Problem]This invention persons develop a new variations numerical orientation 
method, using abinitio molecular dynamics, (i) A superposition phenomenon of fluctuation of each 
process element and a gap of each element was also diagnosed, and a simulation system which also 
enables a check of whether to perform as a sequence of a (ii) request and discovery of a (iii) change 
factor was developed, the further feature of this system — change of a (vi) process — mathematical, a 
point it enabled it to deal with statistically and efficiently, and (v) — a point which added a new 
inference engine, and (vi) — it is in a point which added a new opposite direction system. 
[0023]In a simulation method of physical development in a manufacturing installation of a manufacturing 
process of electronic parts this invention (Claim 1), A manufacturing process factor about the 
characteristic of an electronic industry material formed in said manufacturing installation is made into 
an input value, Are prediction data of the characteristic of an electronic industry material which 
performed a three-dimensional simulation based on this input value, and was formed in a manufacturing 
installation the method of obtaining, and said simulation, It consists of a process performed based on 
molecular dynamics, and a process performed based on a fluid model, At least one side of three- 
dimensional physical quantity outputted as a simulation result performed based on said molecular 
dynamics and three-dimensional physical quantity outputted as a simulation result performed based on 
said fluid model, A simulation method of a manufacturing process of electronic parts which being 
transmitted to another side and using for a simulation of a there is provided. 

[0024]In above-mentioned Claim 1, this invention (Claim 2) from said three-dimensional molecular 
dynamics simulation to a process simulation of said three-dimensional fluid model. An electronic-parts 
manufacture process simulation method, wherein physical quantity, such as a viscoelasticity coefficient, 
a stress constant, a diffusing constant, and Cij, is transmitted is provided. 

[0025]In above-mentioned Claim 1, this invention (Claim 3) from a process simulation of said three- 
dimensional fluid model to said three-dimensional molecular dynamics simulation. An electronic-parts 
manufacture process simulation method, wherein physical quantity, such as a viscoelasticity coefficient, 
a stress constant, a diffusing constant, and Cij, is transmitted is provided. 

[0026]This invention (Claim 4) provides an electronic-parts manufacture process simulation method, 
wherein said electronic parts are semiconductor devices in above-mentioned Claim 1. 
[0027]In above-mentioned Claim 1, this invention (Claim 5) is further provided with a stage of inputting 
fluctuation of said preset value, and provides an electronic-parts manufacture process simulation 
method, wherein said molecular dynamics simulation performs variations calculation based on said input 



value and its fluctuation. 

[0028]This invention (Claim 6) provides an electronic-parts manufacture process simulation method, 
wherein said molecular dynamics simulation is performed in combination of a molecular dynamics 
method and a density functional method in above-mentioned Claim 1 . 

[0029]This invention (Claim 7) provides an electronic-parts manufacture process simulation method, 
wherein said molecular dynamics method is based on the Takashina round partial differential method in 
above-mentioned Claim 6. 

[0030]This invention (Claim 8) provides an electronic-parts manufacture process simulation method, 
wherein said Takashina round partial differential method is dealing with a three-body problem in above- 
mentioned Claim 7. 

[0031]This invention (Claim 9) provides an electronic-parts manufacture process simulation method, 
wherein said molecular dynamics simulation is performed using experiential potential in above- 
mentioned Claim 1 . 

[0032]This invention (Claim 10) provides an electronic-parts manufacture process simulation method 
characterized by a thing in which said experiential potential is ATTA potential, and to do in above- 
mentioned Claim 9. 

[0033]This invention (Claim 11) provides an electronic-parts manufacture process simulation method 
characterized by a thing in which said experiential potential is BKS potential, and to do in above- 
mentioned Claim 9. 

[0034]Namely, according to above-mentioned Claim 1 thru/or the electronic-parts manufacture process 
simulation method according to claim 1 1. A motion of an atom is anew checked by verifying a finer 
portion with a process simulation performed [ calculation / of process induction stress of an extensive 
field of several micron level ] in a result based on molecular dynamics by a process simulation 
performed based on a fluid model. More more exact information can be acquired by performing a 
process simulation of a fluid model again and repeating this if needed. 

[0035]A process of this invention (Claim 12) measuring the characteristic of an electronic industry 
material formed in a manufacturing installation in a manufacturing process of electronic parts, and 
obtaining real observational data, A process of obtaining prediction data of the characteristic of an 
electronic industry material which made an input value a manufacturing process factor about the 
characteristic of an electronic industry material formed in said manufacturing installation, performed a 
three-dimensional simulation based on this input value, and was formed in said manufacturing 
installation, By a process of carrying out comparison assay of said prediction data and the real 
observational data, and said comparison assay. When a significant difference is accepted between a 
preset value of said manufacturing process factor, and said manufacturing process factor guessed from 
said real observational data, Have a process of carrying out correction processing of said manufacturing 
process factor, and said simulation, It consists of a process performed based on molecular dynamics, 
and a process performed based on a fluid model, At least one side of three-dimensional physical 
quantity outputted as a simulation result performed based on said molecular dynamics and three- 
dimensional physical quantity outputted as a simulation result performed based on said fluid model, A 
manufacturing method of electronic parts of electronic parts which being transmitted to another side 
and using for a simulation of a there is provided. 

[0036]A significant difference between said manufacturing process factors it is guessed from said real 
observational data in above-mentioned Claim 12 that this inventions (Claim 13) are has crossed 
tolerance level, When correction processing is impossible, an electronic-parts manufacture process 
simulation method processing as inferior goods is provided. 

[0037]This invention (Claim 14) provides a manufacturing method of electronic parts currently carrying 
out clustering of said manufacturing installation in above-mentioned Claim 12. 
[0038]This invention (Claim 15) provides a manufacturing method of electronic parts, wherein said 
manufacturing installation which carried out clustering is single wafer processing in above-mentioned 
Claim 14. 

[0039]In above-mentioned Claim 12, this invention (Claim 16) from said three-dimensional molecular 
dynamics simulation to a process simulation of said three-dimensional fluid model. An electronic-parts 
manufacture process simulation method, wherein physical quantity, such as a viscoelasticity coefficient, 
a stress constant, a diffusing constant, and Cij, is transmitted is provided. 



[0040]In above-mentioned Claim 12, this invention (Claim 17) from a process simulation of said three- 
dimensional fluid model to said three-dimensional molecular dynamics simulation. An electronic-parts 
manufacture process simulation method, wherein physical quantity, such as a viscoelasticity coefficient, 
a stress constant a diffusing constant, and Cij, is transmitted is provided. 

[0041]This invention (Claim 18) provides a manufacturing method of electronic parts, wherein said 
characteristics are a film formed on a semiconductor substrate, a semiconductor substrate surface, and 
the characteristic about at least one inside a semiconductor substrate in above-mentioned Claim 12. 
[0042]This invention (Claim 19) provides a manufacturing method of electronic parts, wherein 
comparison assay with said prediction data and real observational data is performed in real time one by 
one in above-mentioned Claim 12. 

[0043]This invention (Claim 20) provides a manufacturing method of electronic parts, wherein a process 
of carrying out correction processing of said manufacturing process factor is performed in real time one 
by one in above-mentioned Claim 12. 

[0044]This invention (Claim 21) provides a manufacturing method of electronic parts, wherein said 
molecular dynamics simulation is performed in combination of a molecular dynamics method and a 
density functional method in above-mentioned Claim 12. 

[0045]This invention (Claim 22) provides a manufacturing method of electronic parts, wherein said 
molecular dynamics method is based on the Takashina round partial differential method in above- 
mentioned Claim 21. 

[0046]This invention (Claim 23) provides a manufacturing method of electronic parts, wherein said 
Takashina round partial differential method is dealing with a three-body problem in above-mentioned 
Claim 22. 

[0047]This invention (Claim 24) provides a manufacturing method of electronic parts, wherein said 
molecular dynamics simulation is performed using experiential potential in above-mentioned Claim 12. 
[0048]This invention (Claim 25) provides a manufacturing method of electronic parts characterized by a 
thing in which said experiential potential is ATTA potential, and to do in above-mentioned Claim 24. 
[0049]This invention (Claim 26) provides a manufacturing method of electronic parts characterized by a 
thing in which said experiential potential is BKS potential, and to do in above-mentioned Claim 24. 
[0050]This invention (Claim 27) provides a manufacturing method of electronic parts, wherein a process 
simulation of a fluid model is performed in above-mentioned Claim 12 using at least one prediction data 
of two or more of said characteristics obtained in said molecular dynamics simulation. 
[0051]This invention (Claim 28) provides a manufacturing method of electronic parts, wherein both said 
molecular dynamics simulation and a process simulation of a fluid model are three-dimensional 
simulations in above-mentioned Claim 27. 

[0052]In above-mentioned Claim 27, this invention (Claim 29) from said three-dimensional molecular 
dynamics simulation to a process simulation of said three-dimensional fluid model. A manufacturing 
method of electronic parts, wherein physical quantity, such as a viscoelasticity coefficient, a stress 
constant, a diffusing constant, and Cij, is transmitted is provided. 

[0053]Namely, according to the manufacturing method of above-mentioned Claim 12 thru/or the 
electronic parts according to claim 29. Calculation of process induction stress of an extensive field of 
several micron level with a process simulation performed based on a fluid model a result, A motion of an 
atom is checked anew, calculating a correction amount of a physical property value by verifying a finer 
portion with a process simulation performed based on molecular dynamics. Therefore, a check of 
whether to perform as a desired sequence and discovery of a change factor are also attained, 
performing diagnosis in real time. 

[0054]A measuring device with which this invention (Claim 30) measures a manufacturing installation of 
electronic parts, and the characteristic of an electronic industry material of said electronic parts of 
having been formed in said manufacturing installation, A computer system which makes an input value a 
preset value of a manufacturing process factor of said electronic parts, performs a simulation based on 
this input value, and obtains prediction data of said characteristic, Have a means to connect said 
computer system with said measuring device and said manufacturing installation, and said computer 
system, Perform said prediction data and real observational data, and comparison assay by said 
comparison assay. When a significant difference is accepted between a preset value of said 
manufacturing process factor, and said manufacturing process factor guessed from said real 



observational data, a manufacturing system of electronic parts carrying out correction processing of the 
manufacturing process factor via a control means of said manufacturing installation is provided. 
[0055]This invention (Claim 31) provides a manufacturing system of electronic parts currently carrying 
out clustering of said manufacturing installation in above-mentioned Claim 30. 
[0056]This invention (Claim 32) provides a manufacturing system of electronic parts, wherein said 
manufacturing installation which carried out clustering is single wafer processing in above-mentioned 
Claim 31. 

[0057]This invention (Claim 33) provides a manufacturing system of electronic parts, wherein said 
simulation is performed based on a molecular dynamics simulation by combination of a molecular 
dynamics method and a density functional method in above-mentioned Claim 30. 
[0058]This invention (Claim 34) provides a manufacturing system of electronic parts, wherein said 
characteristics are a film formed on a semiconductor substrate, a semiconductor substrate surface, and 
the characteristic about at least one inside a semiconductor substrate in above-mentioned Claim 30. 
[0059]This invention (Claim 35) provides a manufacturing system of electronic parts, wherein 
comparison assay with said prediction data and real observational data is performed in real time one by 
one in above-mentioned Claim 30. 

[0060]This invention (Claim 36) provides a manufacturing system of electronic parts, wherein a process 
of carrying out correction processing of said manufacturing process factor is performed in real time one 
by one in above-mentioned Claim 30. 

[0061 ]In above-mentioned Claim 30, further, fluctuation of said preset value is inputted as said input 
value, and this invention (Claim 37) provides a manufacturing system of electronic parts, wherein said 
molecular dynamics simulation performs variations calculation based on said input value and its 
fluctuation. 

[0062]This invention (Claim 38) provides a manufacturing system of electronic parts, wherein said 
molecular dynamics method is based on the Takashina round partial differential method in above- 
mentioned Claim 30. 

[0063]This invention (Claim 39) provides a manufacturing system of electronic parts, wherein said 
Takashina round partial differential method is dealing with a three-body problem in above-mentioned 
Claim 38. 

[0064]This invention (Claim 40) provides a manufacturing system of electronic parts, wherein said 
molecular dynamics simulation uses together a method which used experiential potential in above- 
mentioned Claim 30. 

[0065]This invention (Claim 41) provides a manufacturing system of electronic parts characterized by a 
thing in which said experiential potential is ATTA potential, and to do in above-mentioned Claim 40. 
[0066]This invention (Claim 42) provides a manufacturing system of electronic parts characterized by a 
thing in which said experiential potential is BKS potential, and to do in above-mentioned Claim 40. 
[0067]This invention (Claim 43) provides a manufacturing system of electronic parts, wherein a process 
simulation of a fluid model is performed in above-mentioned Claim 30 using at least one prediction data 
of said characteristic obtained in said molecular dynamics simulation. 

[0068]This invention (Claim 44) provides a manufacturing system of electronic parts, wherein both said 
molecular dynamics simulation and a process simulation of a fluid model are three-dimensional 
simulations in above-mentioned Claim 43. 

[0069]In above-mentioned Claim 43, this invention (Claim 45) from said three-dimensional molecular 
dynamics simulation to a process simulation of said three-dimensional fluid model. A manufacturing 
system of electronic parts, wherein physical quantity, such as a viscoelasticity coefficient, a stress 
constant, a diffusing constant, and Cij, is transmitted is provided. 

[0070]That is, according to the manufacturing system of above-mentioned Claim 30 thru/or the 
electronic parts according to claim 45, it becomes possible to predict various physical problems in real 
time as a diagnostic tool in real time by a simulation by a computer. 

[007 1]A procedure of inputting a preset value of a manufacturing process factor about the 
characteristic of an electronic industry material which this invention (Claim 46) is a program which 
performs a simulation of physical development in a manufacturing installation of a manufacturing 
process of electronic parts, and was formed in said manufacturing installation, A molecular dynamics 
simulation by combination of a molecular dynamics method and a density functional method is 



performed based on this input value, and a recording medium which recorded a program which performs 
a procedure of obtaining prediction data of said characteristic and in which computer reading is possible 
is provided. 

[0072]This invention (Claim 47) provides a recording medium, wherein said molecular dynamics 
simulation is based on combination of a molecular dynamics method and a density functional method in 
above-mentioned Claim 46. 

[0073]This invention (Claim 48) provides a recording medium, wherein said molecular dynamics method 
is based on the Takashina round partial differential method in above-mentioned Claim 47. 
[0074]This invention (Claim 49) provides a recording medium, wherein said Takashina round partial 
differential method is dealing with a three-body problem in above-mentioned Claim 48. 
[0075]This invention (Claim 50) provides a recording medium, wherein said electronic parts are 
semiconductor devices in above-mentioned Claim 47. 

[0076]In above-mentioned Claim 47, this invention (Claim 51) is further provided with a procedure of 
inputting fluctuation of said preset value, and provides a recording medium, wherein said molecular 
dynamics simulation performs variations calculation based on said input value and its fluctuation. 
[0077]This invention (Claim 52) provides a recording medium, wherein said molecular dynamics 
simulation uses together a method which used experiential potential in above-mentioned Claim 47. 
[0078]This invention (Claim 53) provides a recording medium characterized by a thing in which said 
experiential potential is ATTA potential, and to do in above-mentioned Claim 52. 
[0079]This invention (Claim 54) provides a recording medium characterized by a thing in which said 
experiential potential is BKS potential, and to do in above-mentioned Claim 52. 

[0080]This invention (Claim 55) provides a recording medium, wherein a process simulation of a fluid 
model is performed in above-mentioned Claim 47 using prediction data of said characteristic obtained in 
said molecular dynamics simulation. 

[0081]This invention (Claim 56) provides a recording medium, wherein both said molecular dynamics 
simulation and a process simulation of a fluid model are three-dimensional simulations in above- 
mentioned Claim 55. 

[0082]This invention (Claim 57) provides a recording medium, wherein physical quantity, such as a 
viscoelasticity coefficient, a stress constant, a diffusing constant, and Cij, is transmitted to a process 
simulation of said three-dimensional fluid model from said three-dimensional molecular dynamics 
simulation in above-mentioned Claim 55. 

[0083]That is, without the ability of actually working a process unit to also write special preparation and 
expense according to above-mentioned Claim 46 thru/or the recording medium according to claim 57, a 
virtual semiconductor manufacturing system can be built and process-simulation ****** can be 
performed to an inside of a computer. For example, it becomes possible about looking for a new material 
easily and effectively moreover by conducting a thought experiment using various molecular structure 
models, and predicting the characteristic of a device using it. 

[0084]This invention (Claim 9) by a molecular dynamics simulator which gave an abinitio molecular 
dynamics process simulator or experiential potential. . It can set to at least one of two or more of the 
processes of a manufacturing process of a semiconductor device. At least one preset value of two or 
more manufacturing process factors and fluctuation of this preset value are made into an input value, A 
film which performed variations calculation based on this input value, and was formed in a manufacturing 
installation in at least one of said two or more of the processes, A means to obtain a film formed on a 
semiconductor substrate, a semiconductor substrate surface, and at least one prediction data of two or 
more characteristics related with at least one inside a semiconductor substrate, Comparison assay of 
the real observational data produced by measuring at least one of said two or more of the 
characteristics in at least one of this prediction data and said two or more of the processes is carried 
out in real time one by one, A simulator which diagnoses and judges a manufacturing process factor 
possessing a means to perform diagnosis and a judgment of said manufacturing process factor of a 
semiconductor device is provided. 

[0085]A very problem that this invention persons described above about a molecular dynamics simulator 
on the other hand that calculation is slow, And an execution area of the target substance was able to 
introduce and conquer a respectively new idea to a problem that it is very small, was able to make a 
molecular dynamics simulator and a process variation calculation program have been able to unite by a 



new method, and, moreover, was able to pull out the new function which is not until now with this union. 
It explains in order below. 

[0086]It touches somewhat in detail about a problem of calculation speed of (1) first In molecular 
dynamics, so to speak, according to Newtonian mechanics, the equation of motion is stood, this is allied 
over 100 pieces or all 1000 atoms to which its attention is paid, and it solves about an atomic piece 
piece. 

[0087]Calculation of power which becomes a basis of this equation of motion is the most complicated. 
This is the late reason by the conventional calculation technique. Each atom has received power from 
an atom around it. Power is got also from a long distance atom also from a nearby atom. It is necessary 
to take this into consideration. It asks for potential energy Vy generally drawn according to a formula of 

TERESOFU, and based on this, power is computed and it goes. However, the technique of making a 
table for every direction from the former in general about power in which potential is beforehand added 
to each atom based on distance, a direction, etc. was taken. About 100 pieces or the 1000 above- 
mentioned atoms, with reference to a table, a potential size is seen and it goes by this method each 
time. This takes time dramatically. Power cannot be computed and inclination to each direction of a 
potential size changes from potential [ itself ] to power of each direction (refer to drawing 2). Because 
of this, it turns out that it takes time of a computer dramatically. 

[0088]It is because the formula showing potential shape is dramatically complicated as to why it makes 
a number table to this appearance. As stated also above, an angle (theta) with a mutual position of a 
surrounding atom, distance (r) of potential (V), etc. are entangled intricately. It is dramatically difficult to 
differentiate this, therefore this mathematically and to make a function. For example, although theta is a 
function of r, a surroundings time grows into a function of r again. 

[0089]Generally it was dramatically complicated and this invention persons created the new Takashina 
round partial differential type conventionally in a portion by which partial differential was not tried about 
such a case so that an argument of an expression potential as mentioned above may go round, thereby, 
although modification of the formula itself becomes apparently very complicated, calculation is burned 
very much — it came to be able to do It is necessary to have ceased to refer to a table like before. 
[0090]The first concept of this Takashina round partial differential type is shown below. Generally this is 
considered supposing three atoms. Namely, if three atoms are generally described, the rest extends this 
technique and should just go. If there are three atoms, distance with each atom, an angle between the 
atom to which its attention is paid, etc. are computable. By these, potential can be searched for 
uniquely and potential "inclination" can be easily computed in any direction based on this. This 
invention persons were making a partial differential type altogether about a variable patrolled and used, 
allying these, and solving, were right and checked that a differentiated type was made. 
[0091 ]An outline of a system by this invention is explained concretely below. First, the actual condition 
of application of this invention is explained with reference to drawing 1. That is, process unit V is some 
semiconductor manufacturing systems which are actually working. Here, an absorption-of-light 
spectrum is seen by spectrum meter VI. 

[0092]An absorption-of-light spectrum acquired by spectrum meter VI is once stored as the file 4, and 
is sent to the scheduler VII. A predicted value by this invention mentioned later is also inputted into the 
scheduler VII, and these are compared with it by the contrast diagnostic assay part IV. 
[0093]If a comparison result in the contrast diagnostic assay part IV is over an acceptable value, it will 
be removed from a line as what does not fill spec. If it is filling, a part for the change will be sent to the 
process variation calculation unit I as an inference engine. 

[0094]In the process variation calculation unit 1, broad view character which is not fundamentally taken 
into consideration to an atom level is calculated, and it stores in the file 2. As for the contents of the 
file 2, a position vector between atoms and a bond angle are computed by the molecular dynamics 
simulator II. 

[0095]It is inputted into the membrane-characteristics calculation unit III, basic membrane 
characteristics are computed here, and a position vector between atoms and a bond angle which were 
computed by the molecular dynamics simulator II are returned to the file 3. Basic membrane 
characteristics computed in the membrane-characteristics calculation unit III are returned to the 
process variation calculation unit I, and let them be a part of a parameter. 



[0096]Thus, a manufacturing method of a semiconductor device which makes it possible to go on by this 
invention while a request corrects a semiconductor device manufacturing process without a test piece 
as a process is realized. 

[0097]The contents of the molecular dynamics simulator II are important especially actually. As shown 
in a block diagram of drawing 7, specifically, it consists of two portions. The 1st portion Ha uses the 
Takashina round partial differential method introduced by this artificer, and the 2nd portion lib conducts 
analysis in consideration of a high order paragraph more than Miyoshi. Details of analysis in 
consideration of the Takashina round partial differential method or a high order paragraph more than 
Miyoshi are explained one by one henceforth. 

[0098]That is, this invention persons developed art which also predicts various spectra by a new 
computational chemistry theory while they left and returned to pure mathematics theory and 
established a new variations numerical orientation method. Hereafter, the new abinitio molecular motion 
study in this invention is explained. 

[0099]An entire configuration of the new simulation which this invention person etc. developed is shown 
in drawing 3. The new technique of having compounded a molecular dynamics portion (MD) and a 
density functional (DFT) was proposed. Proper use with a molecular dynamics method and a density 
functional method by this invention persons is shown below first, using drawing 3. 

[01 00] As shown in drawing 3, a molecular dynamics (MD) portion occupied a lower half of a flow chart, 
and a density functional portion (DFT) occupies an upper half, a molecular dynamics method does not 
come out OK — the so-called a little big system in limited temperature is treated, and potential 
beforehand searched for from the ultimate cause is used for potential between atoms (ion). On the 
other hand, a density functional (DFT) was made to use, when searching for a ground state of a small 
system. And naturally a portion of DFT set temperature to OK. It asked for a potential portion from a 
ground state of an electron system. A wave function (KS equation) was specifically solved and a partial 
density functional was used together. 

[0101]In a DFT portion, it places in quest of total energy of a ground state of an electron system 
without degeneracy. Total energy is functional [ of electron density (** which does not go back to wave 
function) n (r) ] E [n (r)]. 

[0102]That is, E [n (r)] is expressed by following formula (1201). 

[0103] 

[Equation 1] 

Vz of the 1st paragraph of the right-hand side (r) assumes outer fields, such as a nuclear Coulomb field, 
the 2nd paragraph assumes the Coulomb interaction energy between electrons, and the 3rd paragraph 
assumes electronic one-body approximation, and it is the kinetic energy Ts [n (r)] in case there is no 
interaction between electrons, and is expressed by the following formula (1202). 
[0104] 
[Equation 2] 

The 4th paragraph is the exchange correlation energy Exc [n (r)], and the form is based on partial 
density approximation. The electron density n (r) is expressed by the following formula (1203). 
[0105] 
[Equation 3] 

A ground state solves the KS equation (Kohn-sham equation) showing in the following formula (21) and 

(22) which obtains the variations of E [n (r)] of an equation (1201) by setting it with 0, and it is used for 

it N pieces from the one where the characteristic value E is smaller, since n (r) is contained in the 

formula (121 1) — n (r), psi (r), and effective potential Veff (r) must be solved self-contradicting [ no ] 

(self-consistant). 

[0106] 

[Equation 4] 

["** 2+ V eff (r)] psi(r) =Epsi (r) — (121 1) 

V eff (r) =V 2 (r)+integral — 2n(r')/(|r-r'|) and dr'+deltaE xc [n(r)]/deltan (r) 
— (1212) 

The exchange correlation energy Exc [n (r)] obtained a following formula (1213) by applying and 
integrating partial density approximation (LDA) with the electron density n to energy epsilonxc per 



electron (n). 
[0107] 
[Equation 5] 

E xc [n(r)] =integralepsilon xc (n (r)) n (r) dr — (1213) 

The form of epsilonxc (n) is proposed variously the formula besides Wigner etc. An inner shell electron 
is Vz of a formula (121 1), when treating only a valence electron, without thinking. A nuclear Coulomb 
field is not used as it is as (r), but it sets by the norm preservation pseudopotential which erased the 
singularity in the nuclear position, and is a frog. The number of the ingredient in the Fourier expansion 
was able to be pressed down by accustoming in this way. 

[0108]As total energy, this invention persons used the thing also containing the coulomb potential 
between atoms (ion) in order to expect ****. That is, it is a following formula (1221). 
[Equation 6] 

{Rj} is coordinates of the core of ion and {alphanu} is external constraints (volume omega, the strain 

epsilonLnu, etc.) here. The 1st paragraph of the right-hand side of a formula (1221) is Ts of the formula 
(1201) of DFT to the right-hand side. It is [n (r)] and the remainder of the formula (1201) right-hand 
side is one copy of U of a formula (1221). It is expressed with V eff (r) from the formula (121 1) of one- 
body approximation like the case of DFT, and the portion is Vz too. It uses by pseudopotential, replacing 
(r). 

[0109]A new view was put in here. That is, it considers that E is potential energy and the virtual kinetic 
energy shown in a following formula (1222) is introduced into others. 
[0110] 
[Equation 7] 

Therefore, a Lagrangian is expressed by the following formula (1223). 
[0111] 

L=K-E — (1223) 

Although Mj in a formula (1222) was the mass of the true atom (ion), mu and munu is virtual mass and 

was changed according to the purpose like the after-mentioned. The regular rectangular cross 
conditions shown in a following formula (1224) as a binding condition [Equation 8] 

integralpsi*(r, t) psi R (r, t) dr=delta jk — (1224) 

Since it ****, a variations is taken, although undetermined-multipliers lambda jk was introduced and a 

following formula (1231) was added to L, when making the Eulerian equation. As a result, a following 
formula (1232), (1233), and (1234) are obtained. 
[0112] 
[Equation 9] 

The temperature T of the kinetic energy of a formula (1222) corresponds. If it says by the relation of 
sigma1/2 MjRj2 because the degree of means has temperature each one by the energy equipartition law 

of classical statistical mechanics, it will be a temperature that it is not virtual and actual, psii and 
alphanu can be controlled or reinforced by the size of mu and munu. 

[01 13]For example, if it is referred to as T->0 from an elevated temperature as mu«M I? psii will change 

with Rj stopped and the ground state of an electron system will be acquired. The left side of an 

equation (1231) is set to 0 at this time, {lambda ik } which makes the linear combination of psii suitably by 

the equation (1221) and the cautions under it since it becomes a following formula (30) is diagonal-line- 

ized, and the right-hand side obtains a KS equation (1235). 

[0114] 

[Equation 10] 

[** 2 -V eff (r)] psijGO+sigmaps^ (r) — (1235) 

That is, it means obtaining the same result as having solved the KS equation by DFT by simulated 
annealing. However, although the state of the temperature T is made, since it was not based on the 
Monte Carlo method but movement of the virtual dynamical system is followed, it is called 
dynamicalsimulated annealing (DSA). 



[01 15]When integrating with a formula (1231), psii cannot necessarily be treated as a variable as it is, 
but the expansion coefficient (that is, Fourier coefficient) which developed it by the plane wave can be 
treated as a variable. In order not to increase the number of this expansion coefficient recklessly, it is 
required like the case of DFT to erase the singularity in a nuclear position at pseudopotential. When the 
periodicity of a crystal can be used especially, the linear combination of psii which diagonaHine-izes 
above (lambda jk } is described as phink. k is a wave number vector in a Brillouin zone, and n is a zone 

index. 

[0116]phink should fill a following formula (1241). 
[0117] 

[Equation 1 1] 

muS nk (r, t) =[** 2 -V eff (r)] phi nk (r) +lambda nk phi nk (r, t) — (1241) 

lambdank is a characteristic value of a procession {lambda ik } and is an energy level, phink (r, t) is 

developed and a following formula (1242) is obtained. 
[0118] 

[Equation 12] 

If pseudopotential is local (this condition generally is not fulfilled), effective potential V eff will also be 

developed and a following formula (1243) will be obtained. 
[0119] 

[Equation 13] 

A formula (1242) and a formula (1243) are substituted for a formula (1241), and a following formula 

(1251) is obtained. 

[0120] 

[Equation 14] 

It is sigmaexp if it sets with K=K'+K" in the time of the last of the right-hand side. Since [i(K+G) r] 
becomes common to all paragraphs, it removes, and a following formula (1252) is obtained. 
[0121] 

[Equation 15] 

muT n. K +G to =[-|K + G| 2 + lambda nk ] C n and K+Q (t) -sigmaV G _ G ' (t) C n and K+Q ' (t) =-[|K+G| 2 +V 0 (t)- 
lambda J C n and K+G " si g™ V G-G « C n and K+Q ' (t) - (1252) 
Frequency is defined here, as shown in a following formula (1 253), [Equation 16] 
omega=(|K+G| 2 +V 0 (tHambda nx ) (/mu) 1/2 — (1253) 

If the left side of a formula (1251) is changed to difference, a following formula (1261) will be obtained. 
[0122] 

[Equation 17] 

C n K+Q (t+deltat) =2cos(omegadeltat) C p K+G (t)-C n K+G (t-deltat) 
-2 — [1-cos (omegadeltat)] {sigma G 'V G _ G '(t) and C n K+Q ' (t)} 
— (1261) 

Since this formula is clearing up the integration of the vibrating part analytically, it can decrease 
computational complexity by making deltat larger than a numerical method purely. 
[0123]The Coulomb potential is as follows. That is, this invention person found out being divided into 
the 4th paragraph, when the Coulomb potential was decomposed. Although there was only the 3rd 
paragraph in particular conventionally, it checked that there was the 4th paragraph newly. These 
methods are explained one by one. 

[0124]This invention persons began to solve as a start of a strict formula first from the basic equation 
which considered the dielectric constant. It began to solve from the following formula (1262) which is a 
basic equation first. 
[0125] 

[Equation 18] 

a conductor — epsilon=infinity and in the case of the vacuum epsilon= 1, the Coulomb potential differs, 
L is taken by 1 ** of a unit crystal (cube), sigma is taken within a unit crystal, and Zi and ri are 



electrifications and the positions of the i-th particle, respectively. This is because polarization arises in 
the inner surface of a ball by electrification in a ball. Although the layer of electric doublet is made in 
the inner surface of the ball which is not a conductor, the last paragraph of the above-mentioned 
formula serves to negate the effect exactly. Since the method of Ewald gives the thing of the left side 
and epsilon=infinity, it must add the last paragraph of the above-mentioned formula in quest of the 
value within a vacuum. Only a result is hung up. 

[01 26]A following formula (1271) shall express coulomb potential. [Equation 19] 

In the formula (1271), the atomic number in a unit crystal, and atom [ 1st ] position and electrification of 
N in it were ri and Zi, and n is a vector which specifies a unit crystal and the thing which shifted it 
periodically, and as shown in a following formula (1272), it was set up. 
[0127] 

n=n x Z x +n yV n z Z z -(1272) 

It is expressed. Zx, Zy, and Zz are the vectors of Tooru of x, y, and the direction of z here, respectively, 
and nx, ny, and nz are the integers ranging from -infinity to -^infinity, respectively (in the case of a bulk 
crystal). It shall mean that ' of sigma' removes j=i at the time of n= 0. 

[0128]This invention persons introduce here new F function shown in a following formula (1273). 
[0129] 

[Equation 20] 

In the above-mentioned formula (1273), G (r, t) is a periodic function oft. This found out that it could 
express by Fourier series. 

[0130]G (r, t) can change, as shown further below (1281). 
[Equation 21] 

In the above-mentioned formula (1281), m is a reciprocal lattice vector and is expressed by the 

following formula (1282). 

[0131] 

[Equation 22] 

m=m x -(1/L x ,0,0) +m y - (0, 1/L y , 0) +m z - (0,0,1 /L z ) — (1282) 

The index mx, my, and mz are the integers ranging from -infinity to infinity, respectively. At the time of 
m= 0, it is exp. [ — ] Since it is =1, it is sigmaZi =0, and if the total electrification in a unit crystal is 0, a 
paragraph of m= 0 will disappear. If a paragraph of m= 0 is removed beforehand, it will become a 
following formula (1283). 
[0132] 

[Equation 23] 

This invention persons divide an integral range by k, use two forms of G (r, t) properly, and get a 

following formula (1291). 

[0133] 

[Equation 24] 

This invention persons used the following formula (1292) as formal here. 
[0134] 

[Equation 25] 

Thereby, the Coulomb potential to begin is expressed by the following formula (1293). 
[0135] 

[Equation 26] 

' seal at the upper right of sigma' means removing j=i at the time of n= 0. 

[01 36]A following formula (1301) is materialized. 

[0137] 

[Equation 27] 

Since this formula (1301) does not contain r, it is not related to power. This formula (1301) can change, 

as further shown in a following formula (1302). 

[0138] 

[Equation 28] 

When the above result is summarized, the Coulomb potential is expressed by following formula (131 1) - 
(1315) after all. 



[0139] 

[Equation 29] 

n is n=nx, +ny, and (0, Ly, 0) +nz- (0, 0, Lz) here, and m is m=mx, +my f and (1 [ 0 and ]/Ly, 0) +mz - (1 
[ 0, 0, and ]/Lz). The point of the above-mentioned expression in good order is that decrease by the 
factor of erfc in phi (1), and the paragraph of sigma declines quickly by the factor of exp by phi (2) to 
the paragraph of sigma of a basis declining only to the order of a reciprocal in the paragraph of sigma. 
Since it is effective for reverse at slowness and fastness with attenuation by phi (3), it is necessary to 
choose suitable kappa which can balance both. These are the results of calculating contribution of 
Coulomb force sequentially from the near place of distance, and are results in case the circumference 
is moreover a conductor. 

[0140]according to a case where the circumference is a vacuum — already — the 1st paragraph is 
added. Especially this is the portion which was not considered conventionally. 

[0141]Generally in molecular dynamics, interatomic potential is the most important quantity. This also 
has fairly many that are simply described only in two bodies under the present circumstances. Time 
requires molecular dynamics calculation dramatically and a simple thing is being used for this for 
potential for the decrease of a clause of time, thus, if potential is made simple, the calculation itself will 
become quick, but it cannot be a thing reflecting the actual condition. 

[0142]Generally important interatomic potential can be described to have described above by three 
bodies. The implications of this are as follows. Namely, it is effective in starting further first atomic (i) 
and an operation with the second atom (j) through the second (j) from the third atom (k). A thing 
illustrating this is drawing 4. Make r into distance between particles and let theta be a particle angle so 
that it may be shown in a figure. A relation of principal and accessory of these amounts of many is as 
being shown in drawing 5. It turns out that it is dramatically complicated so that it may understand 
especially by drawing 5. Although what differentiated potential generally becomes power, calculation 
becomes very complicated so that it may understand now. Here showed the position price of the 
conventional calculation technique and this invention once again using drawing 6. 

[0143]This invention persons develop newly the Takashina round partial differential method in drawing 6, 
and enabled it to calculate it with high precision and at high speed by this. A actual comparative 
example is put into below and it is described in it how quick methods, such as this invention person, are. 
First, it cannot touch about simple potential calculation between two bodies in drawing 6 here. Because, 
since accuracy which almost bears for describing physical development of semiconductor LSI at use 
was not obtained, it decided not to argue here. On the other hand, this invention persons created a 
program of a table reference method as a conventional method by potential which put in a three-body 
problem. According to this, a table was created about ry, thetay k , and r k . Actually, as block data, 

arrangement was used and it stored. And a thing when it had shifted from this reference item was 
considered, and a complement program was also created. The cubic spline method was used for a 
complement. It was made a program which can quote a table about a certain value, i.e., ry, thetay k , and 

r k . However, desired conditions (ry, thetay k , and r k ) need to explore into which procession within 

arrangement it goes. Although this invention persons created a program, distributing mind so that it may 
change as at high speed as possible, they needed a portion which divides too when based on a size of a 
number, and used an if then else sentence there. Thus, while finding out into which frame of a table 
desired conditions (ry, thetay k , and r k ) fit, it is main ****** about a complement in it. Although potential 

in desired conditions (ry, thetay k , and r k ) was able to be found now first, in fact, I want power in this 

point. In a portion which carried out the minute change, this invention persons searched for potential 
also about one more point using same procedure, and asked for differentiation, i.e., power, from both 
difference using an idea of an average rate of change. 

[0144]Very many if then else procedure was used so that it might understand now. In an electronic 
computer, he is comparatively poor at if then else procedure, and it is slow. In fact, when a referred to 
type program which this invention person etc. created was used, it investigated for details only about a 
portion which calculates power. As a result, by a monocrystal substrate which consists of 1600 Si 
atoms, although it was a case with an external pressure of 1 atmosphere in the absolute temperature 
300K, a referred to type program took about 350 times as much time. It turns out that the Takashina 



round partial differential method which this invention person etc. propose is more advantageous so that 
it may understand now. In a referred to type program, whenever atoms differ, it is necessary to create a 
data table each time. 

[0145]Next, a material property value etc. are extracted from molecular dynamics how, and it is also 
related for how putting it into a general-purpose simulator", and explains briefly. 

[0146]Usually, evaluation art of all sorts is used in an LSI process. A property value appraisal method of 
a film by infrared absorption, stress evaluation in a film according to the Raman spectrum further, 
structural analysis of a film by X-rays, etc. occur. There is not only these optical things but electric 
evaluation art, such as a membranous dielectric constant. 

[01 47]It decomposed even into movement of an individual atom theoretically once wholeheartedly, and 
this invention persons pursued a principle of each evaluation technique in detail. Using molecular 
dynamics which this invention person etc. moreover developed, a motion of an individual atom was able 
to be pursued, and each evaluation quantity was able to be computed theoretically for the first time, for 
example, it was able to display directly on a form like a spectrum. Thereby, it is a reason for having 
obtained progress of big art which is not in the former that quantitive comparison with an actual 
measurement can be performed etc., for example. 

[0148]this invention person etc. describe even pursuit of atomic movement by the first theoretical 
means to conquer, and derivation of the concrete amount of measurement evaluation below. Although 
there are some things, fundamentally, the fundamental concept which this invention person drew is as 
follows. That is, a moment of an atom to which its attention is paid is calculated. And this is integrated 
in a window of time to pay one's attention, the Fourier transform is carried out to after an appropriate 
time, and an inherent spectrum is obtained. This is a big main point. 

[0149]Next, Si is shown [ experiential ] for how this invention person changed the potential into power- 
by being related potentially in an example. Although this is an example, if it does not take into 
consideration an effect not only by Si but an ionic bond, the Takashina round partial differential 
described here can be used for it as a three-body problem, for example like [ potential of an oxide film ] 
a problem of this Si. What kind of potential is good examined this invention persons about Si. 
[0150]Potential used here is based on Tersoff. (Phys.Rev.Lett.56.632 (1986), Phys.Rev.B37-6991 
(1988)). 

[0151]Total potential [ according to now and Tersoff] about i-th Si, [Equation 30] 
It can come out and describe. Since three bodies have described this, it is Vy!=Vji in the above- 
mentioned formula (0001) that this invention person etc. want to call attention. When the location 
number of Si atom to which its attention is paid is set to i and other particle numbers of the 
circumference of it are set to j, the above-mentioned Vy is, [Equation 31] 

V U =f o W + Va Org)] ~ (0002) 

It comes out. r is the distance between particles here, fc (ry) is a cutoff function, fR (ry) shows repulsive 

force and fA (ry) shows attraction, ay are a cutoff coefficient in consideration of the coordination 

number, and the cutoff coefficient as which by also considered the coordination number. Three-body 

expression will be performed using a parameter expressing the coordination number. Since a basic form 
is the same in other things, this shows below here, mixing thought of the first round technique for this 
invention person etc. to develop. 

[0152]fR (ry) and fA (ry) are what changed Morse type potential, and are fR(ry) =Aexp (-lambdalr) and 
fA(ry) =-Bexp Hambda2r). 

[0153]Among this, lambdal and Iambda2 are constants, and that size is a reciprocal of a value about 
interatomic distance. 

[0154]By the way, the cutoff function fc (ry), [Equation 32] 

V ii" f r ( r ii) [a ii Aexp(-lambda 1 r i! )-b ;i Bexp (-lambda 9 r ;; )] 
U c U U I u u ^ u 

f c (r) =1 (r<=R-D) 

f c (r) =1/2-1 / 2sin [pi/2(r-R) /D] (R-D<r<R+D) 
f (r) =0 (r>=R+D) 

C 



— (0003) 

It comes out, and it is, and R chooses the size here so that only the first adjacent zone of the target 
structure may usually be included. According to this invention person's detailed examination, the value 
is about 2-3A. Next, although by will show the real coordination number, it uses the above-mentioned 

cutoff function also here. According to Tersoff, the definition is b ij =(1+beta n zeta jj n ) ~ 1/2n — (0004) 
It comes out. Here, [Equation 33] 

zeta jj =sigmaf c (r jk ) g(theta jjk ) exp [lambda 3 3 (r ij -r jk ) 3 ] — (0005) 

It comes out. sigma sign turns by k!=i and j. Since the meaning of zetay is an environmental factor by 
the 3rd atom k entering as it understands here, it is with the case where it sees from i, and the case 
where it sees from j, and sizes differ mutually. That is, it is zetay!=zetajj, and further, as the above- 
mentioned formula (0001) described, it is VJ=Vji. Therefore, it is b j j!=bji. 
[0155]g (theta) is a bond angle factor, [Equation 34] 
g(theta) =1-(c 2 /d 2 )-c 2 /(d 2 +costheta 2 ) — (0006) 

It comes out. Here, theta shall be taken like drawing 4. In asking for theta, it expresses using actual 
rectangular coordinates. namely[Equation 35] 

r.=root {(Xj-Xj) 2 +( yj - yj ) 2 +(z j -z j ) 2 } — (0007) 

It comes out, and it is and r- k is also called for in same procedure. If it does so and an inner product will 
be made into Py, [Equation 36] 

p yk = ( x j- x i) (Xk-XiMyj-Vi) (y k -Vi> + (z k - Zj ) - (ooos) 

It comes out. These are used and it is costhetay k =Py k /(ry-r jk ). — (0009) 

It becomes. Here, the constant in each above-mentioned formula which this invention person etc. used 
is shown. namely[Equation 37] 

R= 3.0 A, D= 0.2 A, A=3264.7ev, B=95.373ev, C= 4.8381, lambda 1 =3.2394A, lambda 2 =1.3258A, 
lambda 3 =lambda 2 , beta= 0.33675, n= 22.956, and d= 2.0417 — (0010) 
It comes out. 

[01 56]this invention person's round portion is explained in detail still more through these preparations. 
From here, power is calculated still more and it goes. It will become power if a potential (0002) formula 
is differentiated from coordinates of a position. namely[Equation 38] 

It is x ingredient of the vector of power which works to ********** and the particles i and j. However, 
it is not so easy to actually ask for it. Below has the device and invention which stood on this invention 
person's etc. mathematical theories, and returned. It is the division called Takashina round partial 
differential. 

[0157](001 1) What is necessary is just to carry out partial differential of the formula described by the 
formula and (0012) the formula concretely. However, it turns out immediately that it is not so simple. 
[0158]namely[Equation 39] 

Modification of the partial differential equation about j and k is as follows similarly. Although the above 
mainly described x, it is necessary to perform it also about y or z. A blank part keeps blank and was 
placed especially here so that correspondence with the above might be known. Such complicated 
modification was not found as long as this invention person etc. investigated. Such complicated 
modification was conquered until now and, for the method reason of ****, it had escaped to the simple 
table reference method conversely. Since such an invention was conquered, this invention persons have 
united molecular dynamics and a general-purpose fluid simulator. 
[0159] 

[Equation 40] 

Each item of these formulas can change as follows. 
[0160] 

[Equation 41] 

**Vy/**r y =**f c . (ry)/**ry and [Aexp(-lambda 1 /r ij )-b jj Bexp(-lambda 2 /r ij )]+f c (r y ) [-lambda 1 Aexp(- 



lambda.! / r i j)+lambda 2 b i jBexp (-lambdag/ry)] 

= Aexp. (-lambda, /r..) -[**f (r..)/**r..-lambda,f (r..)]- Bexp(-lambda,/r ;; ) - [**f (r ii )/**r ii -lambda 0 f GO] 

i U c y ij £. c ij £. \\ c y ij z c ij 

b B — (0016) 

*^ c (r y )/**r y =-pi/4 D-cos[pi/2-(r-R)/D] 

(R-D<r<R+D) 

**f ( r )/** r ..=0 (r>=R+D) 

C IJ IJ 

— (0017) 

**V y /**zeta y =. (**V y /**b y ) -. (**b ij /**zeta jj ). = -Bfc. (r y ) exp- (Iambda2r a ) - 1/2n " 1 beta n zet aij n ~ 1 =Bfc 

(r y ) exp-(lambda2r y ) b^betazeta^ n /(-1/2n) (1+beta n zeta a n ) [2{1+(betazet aij ) n } zeta y ] 

— (0018) 

**zeta ij /**r jj =3lambda 3 3 sigmaf c (r ik ) g(theta yk ) exp(lambda 3 3 (r |:j -r ik ) 3 ) (r y -r jk ) 2 — (0019) 
[Equation 42] 

**zeta y /**r |k =df c . (r y )/dr y and g(theta jjk ) exp (Iambda3 3 (r y -r ik ) 3 )-3lambda 3 3 f> ik ) f c (r y ) (r--^) 2 - 
(0020) 

**zeta ij /**costheta jjk =f c . (r jk ) exp. (Iambda 3 3 .) (rg— r- k ) 3 dg (theta jjk /dcostheta jjk f c (r jk ) exp(lambda 3 3 (r jj - 

r ik ) 3 ) [2c 2 costheta |Jk /{d 2 cos 2 theta ijk } 2 ].) — (0021) 
**r ij /**x i =(x-x j )/r jj =**r jj /**x i — (0022) 
**r ik /**x i =(x j -x k )/r jk =**r jj /**x k — (0023) 
[Equation 43] 

** — cos — theta — yk __ ~ /~ ** — x — __ . __ — = — one — /(r y r jk ) - ** — P — jjk __ — 
/-- ** — x — - _ . __ " + ~ — P — __ jjk __ — (— one — /(r jk ) — ** — /-- ** — x — __ . ~ [— 
one — /(r jk ) — ] — " + ~ — one — /(r y ) ~ ** — /-- ** — x — _ . „ ~ [— one — /(r jk ) — ] — } 

= -(x r x k+ x rXj )- P ij J(x j x j )/(r ik r y 3 Hx r x k )/(r ik r y 3 ). = 1 /(r jk ) - {(x- Xj ) /(r y )- /(x,-^) (r jk ) costheta y } 

— (0024) 

**costhota Bk /** X| =-1/(r B ) and Kxj-x k )/(r jk ) -(x-xpAr^) costhetay} — (0025) 
**costheta jjk /**x k =-1/(r jk ) and {(x- Xj )/(r y ) -(x r x k )/(r ik ) costhetay} — (0026) 

L (r^r./^t, V, **V/**t) =1/2sigmam(**r j /**t) 2 -U(r j )+1/2M(**V/**t) 2 -P E V-(**L(q j , qj)/**q.) d 
(**L/**q\)/dt=0 — (0027) 

Cooperation with the above-mentioned process variation calculation unit I is explained to drawing 7 and 
drawing 8 as the molecular dynamics simulator II. There are some which are transmitted as the first 
input a from a file (22) shown in drawing 8, (32), (42), (52), (62), (72), etc. in the molecular dynamics 
simulator II. the contents — a boundary condition, and area size (namely, breadth of a x direction, 
breadth of a y direction, breadth of the direction of z) and stress (namely, sigmax, sigmay, sigmaz) — it 
is distorted (namely, epsilonx, epsilony, epsilonz) and is temperature etc. There is the input b auxiliary 
as the second thing if needed, and there are computation time of molecular dynamics, specification of 
ensemble, etc. 

[0161]On the other hand, position vector r (t) comes out as the first output. And it is transmitted to a 
file (23) shown in drawing 8, (33), (43), (53), (63), (73), (83), etc. Membrane characteristics, such as 
optical, electric, and crystallographic characterization, are computed, and it is transmitted to a file (3) of 
drawing 8 by ** and the membrane-characteristics calculation program III. 

[0162]An example of an input sequence of wide use 2 and a three-dimensional process simulator ((la) of 
drawing 8) is shown in drawing 9. wide use 2 and a three-dimensional process simulator ((10) of drawing 
8) will accept it as well as drawing 10 — one input sequence is shown. While contents shown in drawing 
9 also specify a cluster tool used, respectively, for example, there are gas volume, a temperature-up 
sequence, a valve sequence, etc. in details of an oxidation process which is, for example in it. 



[0163]On the other hand, an input sequence of drawing 10 is a case where strange material shown by 
drawing 8 is included. In the case of this drawing 10, it is a ferroelectrics film, and it includes a process 
of carrying out the weld slag of the film, for example. Since there is no detailed data of a ferroelectrics 
film at this time, it flies to the molecular dynamics simulator II, and is once calculated there at it. 
[0164]This invention persons perform a molecular dynamics simulation which gave experiential potential, 
(1) A superposition phenomenon of fluctuation of each process element or a gap of each element was 
also diagnosed, and a simulation system which also enables a check of whether to perform as a 
sequence of (2) requests and discovery of (3) change factors was developed. 

[0165]While proposing physical relationship of a single crystal ferroelectrics film and a Si substrate, and 
position relation with a metal electrode, specifically, an indicator was able to be further shown for the 
first time also about what quality of a thing should be prepared in practice as a single crystal 
ferroelectrics film. While reproducing this single crystal ferroelectric membrane structure with the 
original calculation technique using a computer strictly, a characterization function is given, ****** and 
impression conditions are considered to this, and what should be done with physical relationship with a 
substrate single crystal to it, or [ simultaneously ] an indicator of tolerance level of the defect density is 
shown. 

[0166]This invention persons invented a method of carrying out calculation prediction of various 
phenomena in a system which contains an element variously at high speed and with high precision by a 
new computational chemistry theory. By this, the characteristic of various materials used for a 
semiconductor manufacturing process can calculate now with high precision. 

[0167]As a means to predict conventionally structure and the characteristic of a crystal or a molecule 
that an element is contained variously, an abinitio molecular dynamics method and a ultimate cause 
molecular orbital method have been used. An abinitio molecular dynamics method was mainly used for 
structure of a crystal, or investigation of an electronic state, and a molecular orbital method has been 
used for performing calculation of the structure and electronic state for a cluster which cut down the 
usual molecule and some crystals. However, in these methods, in order to solve a schrodinger equation to each 
electron of a system, huge computation time was required. Therefore, it was actually impossible to have 
incorporated these ultimate cause technique in the technique of needing calculation times of several 
100,000 steps like a molecular dynamics method which gave experiential potential. On the other hand in 
the field of molecular dynamics, A.KRappe and W.A.Goddard in a paper of "Charge Equillibration for 
Molecular Dynamics Simulations" in 1991, The new electric charge calculation technique for a 
polyatomic molecule is proposed. This invention persons extend still more greatly the electric charge 
calculation technique in a polyatomic molecule proposed by A.R.Rappe and others, It changed 
substantially so that it could apply to various elements and crystal systems, and a materials design 
simulation procedure applicable also to PbTi03, a ferroelectric called PZT, and a high dielectric was 
built as well as Si and Si0 2 which are generally used as a material of a semiconductor device. First, the 

technique of Rappe and others is described below. 

[0168]Rappe and others considered the energy E of an isolated atom as a function of ionization 
potential and electron affinity first. That is, in a certain atom A, AO shows a neutrality state electrically. 
If it thinks to the secondary paragraph of Taylor's expansion, [Equation 44] 

An atom is obtained by the energy of a neutrality state substituting zero for QA of a formula (3031) 

electrically. 

[0169] 

E A <°> =E A0~ < 3041 ) 

When an atom +-1 -value-ionizes, energy is acquired by substituting one for QA of (3031), and when it 

ionizes to -1 value, energy is acquired by substituting -1. 

[0170] 

[Equation 45] 

The difference of (3042) and (3041) is set to ionization potential IP, and the difference of (3043) and 

(3041) becomes the electron affinity EA. 

[0171] 

[Equation 46] 

(3045) When it considers it to be positive work to take an electron, - is attached to EA of a formula in 



order to mean that adding an electron had done negative work. If the difference of (3044) and (3045) is 

taken, it will become the electronegativity chiAO. 

[0172] 

[Equation 47] 

A following formula will be realized if the sum of (3044) and (3045) is taken. 
[0173] 

[Equation 48] 

Here, the physical meaning of a formula (3052) is considered. The mimetic diagram of the Coulomb 
interaction at the time of the ionization in the case of a hydrogen atom is shown in drawing 12 thru/or 
drawing 14. As shown in a figure, the difference of IP and EA is Coulomb interaction between electrons 

b. Therefore, it is IP-EA=J AA ° when the coulomb repulsive force between 2 electrons in phiA orbit is 

expressed as JAA0. — (3053) 

JAA0 will be called idempotential. Since orbital shape changes by adding an electron, some values of 
JAA calculated by ultimate cause calculation of Hartree-Fock etc. and a formula (3053) differ, but it 
approximates by a formula here (3053). 

[0174]The relation between Idempotential JAA0 and atomic diameter RAO is considered. Considering 
the easy atomic structure model shown in drawing 15, Coulomb interaction JAA0 between electrons is 
as follows. 
[0175] 

[Equation 49] 

Here, JAA0 is eV unit and R is ** unit. (3053) In order to check the validity of a formula, the relation 
between J value and R was investigated and it collected into Table 1. As shown in Table 1, J was 
mostly inversely proportional to atomic sizes, and it turned out that the diameter of the atom called for 
from the formula (3054) is about in agreement with homopolar bond distance. 

[0176](3051) If a formula is used, energy (3031) of an isolated atom can be written with a following 

formula (3053). 

[0177] 

[Equation 50] 

In the case of a polyatomic molecule, it is necessary to also take the electrostatic energy between 
atoms into consideration. Electrostatic energy can be expressed with a following formula. 
[0178] 

[Equation 51] 

Here, JAB is a Coulomb interaction in case there is 1 electron charge at the center of the atoms A and 
B. JAB is dependent on the distance between AB(s). (3055) When a formula is used, all the electrostatic 
energies can be expressed with a formula (3057). 
[0179] 

[Equation 52] 

Chemical potential chil (Q1, QN) of an atomic scale can be expressed with a following formula. 
[0180] 

[Equation 53] 

The simultaneous equations of KA (N-1) are obtained from the conditions (3059) of an electric charge 

balance. 

[0181] 

chi t =chi 2 = =chi N — (3059) 

namely [Equation 54] 

It becomes. To this, N simultaneous equations including the condition (3062) equation of all the electric 

charges are solved. 

[0182] 

[Equation 55] 

(3061) When it displays by a matrix, it can express with the following formulas (3062). 
[0183] 

[Equation 56] 

(3063) Since a maximum minimum is among the electric charges which each ion can take in solving a 



formula, if an electric charge becomes in the range outside, fix the electric charge to a boundary value. 
(3061) It separates into the electric charge which fixed the formula to the boundary value, and the 
electric charge which is not so, correct a formula (3063), dispel a matrix, and make it converge. Then, 
dispersion in an electric charge is calculable. 

[01 84] All that matters in this method is the determination of JAB. When the distance RAB with the 
atoms A and B is large [Equation 57] 

Although it is easy to come out, it is necessary to consider repulsive force in distance with which an 
electron cloud laps. Here, electron density will be approximated by a paragraph of the one Slater orbital. 
That is, in ns, np, or an atom that is nd, an outside valence orbital presupposes that it consists of the ns 
Slater orbital by which the following mold was standardized in any case. 
[0185] 

[Equation 58] 

(3065) If it asks for atomic average size using a formula, [Equation 59] 

It becomes. Valence orbital index zetaA is chosen with the following relation. 

[0186] 

[Equation 60] 

lambda is a parameter for adjustment, and it is contained in order to calculate the difference of the 
average atomic sizes and covalent-radius RA of a crystal which were given by the formula (3066). 
[0187]The scaling parameter lambda is set up using an alkaline metal halide molecule as shown in 
drawing 16. An alkaline metal halide molecule is set to MX, and use M as an alkaline metal (Na, K, Rb, 
Cs), and let X be halide (CI, Br, I). When each electric charge is expressed as QM and QX, it is 
Q tota| =Q M +Q x =0 thereforeQ x =-Q M from a formula (3062). — (3068) 

(3061) From a formula [Equation 61] 

(3069) If a formula (3068) is substituted for a formula, [Equation 62] 

( J MM- J MX- J XM +J XX> Q M =chi ° x" chi ° M " < 3071 > 
By the way, it is J^x^XM c ' ear ' y - — (3072) 
Because it comes out, [Equation 63] 

(3073) The variable of a formula is only lambda. Dipole moment muMX is calculated using a formula from 

calculated Q (3074), and lambda is determined as compared with an experimental value. 

[0188] 

mu MX =4.80324Q M R MX — (3074) 

RMX is an experimental value of bond distance. (3074) As for the formula, Q becomes, R becomes ** 
per e, and mu has become a unit of a debye. 0.4913 is obtained from comparison with the experimental 
value of a dipole moment as an optimum value of lambda. If this lambda is used (3067), the valence 
orbital index zeta can be found in the value of the table of drawing 1 7 from a formula. 
[0189]About hydrogen, Rappe and others is doing different handling from other elements. 
[0190]As compared with Mulliken and a value on an experiment of Pauling others, it is not in agreement 
in electronegativity of hydrogen. On a Mulliken scale, hydrogen is electronegativity from carbon or 
nitrogen to hydrogen being electropositivity from carbon and being electronegativity in a Pauling scale, 
more slightly than boron. Since the orbit of hydrogen which has combined a problem of this 
electronegativity cannot spread like a free hydrogen ion, it is because it becomes smaller than a value in 
case the execution electron affinity EA is an atom. 

[01 91]Therefore, EAH is made into a variable and chiOH of hydrogen and J0HH are redefined. 
[0192] 

[Equation 64] 

It is chi° H =4.528eVJ° HH (0) =13.8904eV by comparison with an experiment to the least square method. 
— (3076) 

Not an experiment but the comparison with the electric charge calculated from the electrostatic 
potential of HF wave function to chi° H =4.7174eVJ° HH (0) =13.4725eV — (3077) 
It is supposed that it may be used. 

[0193]A boundary value about an electric charge is set up. 



[0194]A valence shell is not appropriate to chi A and J AA out of a range corresponding empty or to the 

limit clearly. For example, in Li, C, and O, an electric charge is restricted at a range as shown in the 
following figure. E A (Q) infinity is taken out of this range. 

[01 95]It judges whether it is within the limits of inequality showing an electric charge which dispelled 
and searched for a matrix (3063) in drawing 18 thru/or drawing 20, and in being outside the range, it 
fixes the electric charge to a boundary value. 

[0196]The above is an outline of the technique of Rappe and others. Rappe and others has applied the 
above-mentioned technique to a molecule of an organic system. This invention persons adopted only a 
concept of chemical potential of an atomic scale which Rappe and others proposed, and newly developed the 
individual atomic charge calculation technique in a crystal system. The electric charge calculation 
technique by this invention persons is described below. 

[0197](3057) The second paragraph of the formula right-hand side means the sum of electrostatic 
energy of all the combination. This invention persons extended to a crystal from a molecule, and tried a 
periodic boundary condition with a method of introduction. Then, when calculating electrostatic energy, 
it discovered that it was necessary to take into consideration not only hundreds - thousands of atoms 
that are used for calculation but electrostatic energy received from a countless virtual atom on the 
outside. That is, attenuation of electrostatic energy was small and he has noticed that influence from a 
distant atom is also steamed and is not made. Then, electrostatic energy from a virtual cell was also 
taken into consideration, assuming that a countless cell exists in the outside of a cell used for 
calculation virtually. 

[0198]In this case, as for the electrostatic energy JAB received from an atom in a formula (3057) virtual 
cell, when a size of and calculation cells IS taken into consideration, RAB is usually a big value of 10 A or 
more. Since it becomes, Coulomb energy and ** may obtain. Therefore, JAB=1/RAB, and 1 It treats to 
**. If a virtual cell is expressed with nu, a formula (3051) will be changed as follows. 
[0199] 

[Equation 65] 

Here, nu= 0 means th e calculation cells themselves and nu!=0 means a virtual ceii. (3078) Apply the Ewaid method 
to the 3rd paragraph of a formula. If the Ewaid method is used, the 3rd paragraph can be searched for as 
follows. 
[0200] 

[Equation 66] 

It is a parameter in which G adjusts a reciprocal lattice vector, omega adjusts the volume of calculation 
cells, and kappa adjusts convergence of the Ewaid calculation. The first paragraph of the formula (3082) 
right-hand side does not include the paragraph of nu= 0 and A=B, and the second paragraph does not 
include the paragraph of G= 0. If it asks for chemical potential chil of an atomic scale from these formulas, 
[Equation 67] 

It becomes. N-1 equation is obtained from conditions of an electric charge balance. 
[0201] 

chi 1 =chi 2 =chi 3 = .... =chi N — (3092) 
(3091) Substitute a formula, [Equation 68] 

(3093) N simultaneous equations which can be expressed with an equation will be solved. 
[0202](3093) A following formula will be obtained if the matrix of CQ=-D is made from a formula. 
[0203] 

[Equation 69] 

This invention persons discovered that what is necessary was just to solve a formula, taking a periodic 
boundary condition into consideration, when the crystal and amorphous substance like Si0 2 were 

treated (3094). (3094) To the electric charge Q, since a formula is nonlinear, it brings the i oop for 
completing Q to a calculation program. 

[0204]This invention persons also paid attention to modeling of the Coulomb integral which is in charge 
of real calculation. JAB was calculated as a two-atom Coulomb integral of the siater function over zetaA 
and zetaB. Coulomb integral JAB. It calculated using the method of ** Roothern. Roothern and others 
uses the sign of being in Coulomb integral JAB. (3101) A Coulomb integral can be expressed with a 



formula (3102) using the parameter defined by the formula. 
[0205] 

[Equation 70] 
Here[Equation 71] 
(alpha, alpha', beta, beta') 

= (1+tau a ) alpha(1-tau a ) alpha'(1+tau b ) beta(1-tau b ) beta' — (3103) 
It comes out. Further,[Equation 72] 

Here, the unit systems of an upper type are atomic units, distance uses bohr and energy is using 
Hartree. 

[0206]Although Roothern and others had given the analysis formula of the Coulomb integral to 1s, 2s, 
and 2p, this invention persons calculated to the Coulomb integral over the atomic orbital after 3s, and 
they chose the optimal orbit so that Si and Si0 2 which are generally used with a semiconductor device 

could be treated. 

[0207]The value of JAB influences electric charge calculation greatly. Therefore, setting out of the 
Slater orbital used when computing J, or the valence orbital index zeta is important. Then, this invention 
persons investigated about the value of J in the case of treating Si-0 and Si-Si. 
[0208]Although it is carrying out "assuming if polarization charge is in ns orbit" Rappe and others, a 
point whether it is optimal taking which orbit into consideration to an element variously is unknown. 
Since the principal quantum number n of Si was 3, it compared by performing calculation a case where a 
model is made when an electric charge of Si was in 3 s, and at the time of assuming that it is in 2 s. At 
this time, this invention persons prepared an analysis type of J in n= 3 separately in accordance with a 
method of Roothern. For example, a following formula was used to Si-Si. 
[0209] 

[Equation 73] 

Drawing 21 is a calculation result of Coulomb interaction J between Si-Si. In this case, even if it uses 
any (2s and 3s), it turns out that there is almost no change in the value of J. However, Coulomb 
interaction J between Si-0 differs greatly by the case where 2s and 3s are used. The electric charge of 
O is shown in drawing 22 at 2 s, and the calculation result of JAB at the time of making a model, when 
Si was in 2 s or 3 s was shown. In response to the influence to which 3 s orbits have spread 
dramatically, the value of J will have [ in / for the electric charge of Si / 3 s ] the minimal value at 
nearly r= 1.6A. Since the bond distance of Si-0 was in this hit exactly, J estimated very small and it 
turned out that a result and any electric charge of Si and O will be evaluated too little [ very ]. The 
same might be said of the case of Si-H. Therefore, as a model of the electric charge of Si, if it dares to 
assume by 3 s orbits, the estimate of J will completely lack accuracy. From this, it was only in order to 
take the breadth of an electric charge into consideration to use the Slater orbital, and it was 
understood that it is optimal as a model to use the Slater orbital for 1 s or 2 s. This was the same also 
to the element of the atomic number after Si, and when the orbit it spread after [ whose ] 3s was used 
as the model, it turned out that the calculation precision of Coulomb force when interatomic distance is 
near gets very bad. 

[0210]This invention persons also paid attention on a program to amendment to hydrogen. (3031) From 
a formula, [Equation 74] 

Chemical potential chiHI of an atomic scale when hydrogen amendment is carried out is calculated. 
First, the formula of the energy of the polyatomic molecule corresponding to a formula (3057) is, when it 
is M hydrogen among all the N atomic numbers, [Equation 75] 

Here, when the sign of sigma' of the first paragraph of the right-hand side of a formula (31 14), the 
second paragraph, and the fifth paragraph asks for each sum, it means that H is removing. Chemical 
potential chiHI of an atomic scale, [Equation 76] 

(3121) The fourth paragraph of the formula right-hand side can be unraveled as follows. 
[0211] 

[Equation 77] 

(3122) A formula shows that the sum of the fourth paragraph of (27) types and the fifth paragraph 
means the Coulomb interaction between atoms between all the atoms other than HI. Therefore, 
chemical potential, [Equation 78] 



(3123) Change will be added to a matrix (3063) at the time of hydrogen amendment using a formula. 
Now, a formula (3123) will be substituted for a formula (3059) supposing chi 1 is atoms other than 
hydrogen, [Equation 79] 

As compared with a matrix (3063), it is about JiiO of hydrogen. [Equation 80] 
It is alike and will change. 

[0212]Next, hydrogen amendment of the Coulomb potential JAH between atoms is considered. The siater 
orbital function in the case of hydrogen is a formula (3075), [Equation 81] 

When hydrogen is contained, after searching for an electric charge, a formula (3126) is used, JAH is 
corrected, and an electric charge is searched for again. And this is repeated until it is completed by the 
value of an electric charge. Under the present circumstances, when the line of 1= 1 was hydrogen, this 
invention persons created the algorithm so that it might change for the line which is not hydrogen and a 
matrix might be created. 

[0213]This invention persons also paid attention to the treatment of an electric charge boundary 
condition for the concrete changing method of the matrix. 

[0214](3061) Divide a formula into electric charge QC which is not so, and consider it to be the electric 

charge QB fixed to a boundary value. 

[0215] 

[Equation 82] 

(3062) When a formula is changed into a formula (3131), the element of a matrix is at the ** time, 
[Equation 83] 

A matrix is changed like an upper type and an electric charge is stored within the limits of a boundary 
condition. When the element of 1= 1 was set as a boundary value at this time, it was made to choose the 
line so that the line of a matrix might be replaced and the line to the element by which boundary-value 
setting out is not carried out might be set to 1= 1 (pivoting). 

[0216]They computed the IR spectrum by this invention persons having added this to the molecular 
dynamics method, having calculated the motion of an individual atom one by one using the above- 
mentioned electric charge numerical orientation method newly developed in the crystal system and the 
amorphous system, having calculated the dipole moment every moment, having recorded it, and having 
done the Fourier transform of this. Drawing 23 summarizes the above procedure. 

[0217]On the other hand, in the gate oxide used for DRAM after 16G, or the tunnel oxide film of NAND 
EEPROM, high reliability-ization has been important SUBJECT also compared with the former. In order 
to use an oxide film under the severe condition of maintaining insulation, in EEPROM especially, sending 
tunnel current under a high electric field, it is important to search an atom and an electronic level for 
the limit as an insulating material. 

[0218]From such a situation, a search experiment with an atom level of Si/Si0 2 interface structure or 

Si0 2 network structure, an experiment of a trap mechanism of Si0 2 , prediction by theoretical 

calculation, etc. are advanced in various research institutions. However, like this invention, relation of 
network structure and IR full spectrum was associated with an atom level, and there was no technique 
to be used in the former, an experimental result which associates an asymmetric stretch peak of Si-0 
of an IR spectrum, and stress of an oxidation interface slightly — a report — now, it was only. 
[0219]Only an argument from a vibration analysis of Si-O-Si molecular structure which simplified an 
actual network dramatically is made using a "central force model" also as an interpretation of an 
experimental result, It did not have a discussion in consideration of all distribution of distribution of 
amorphous characteristic Si-0 bond length, Si-O-Si, and the degree of O-Si-O bond angle at all 
conventionally. In this invention, correlation with an electrical property and structure was acquired from 
a different viewpoint from a central force model by taking into consideration these amorphous 
characteristic various distribution. 

[0220]Here, it adds about a point which this invention persons conquered. In using a classic molecular 
dynamics method, this invention persons examined the potential, and also they added improvement 
further. Although this invention persons tried use of TTAM potential developed by Mr. **** and others 
at first, they found out a fault that Si-0 asymmetric stretching vibration was dramatically expressed to 
the low wave number side. Then, BSK potential developed by van Beest etc. was used. As for their 
potential, an elastic peak was expressed with a wave number lower in addition than an actual 



measurement although a peak appeared in the number side of Takanami also in a TTAM potential twist. 
[0221 ]In order that this invention persons might bring this peak value close to survey, they removed 
conditions of electric charge regularity currently used by these existing potential, and incorporated 
distribution of an electric charge when a Si-O-Si angle etc. are situated very much into classic 
molecular dynamics. Calculation of dealing with hundreds of atoms structure at once, and performing 
electric charge distribution calculation as an actual problem, under the present circumstances although 
calculation of this electric charge distribution is ideally computable using a molecular orbital method It is 
actually impossible even if it uses a supercomputer. 

[0222]Then, this invention persons performed electric charge distribution calculation using a technique 
proposed by A.R.Rappe and W.A.Goddard III. Citation papers are "Charge Equillibration for Molecular 
Dynamics Simulation", J.Phys.Chem., 95 and 1991, and p.3358-3363. This invention persons developed 
expression for their paper to reference, after correcting an error in a paper, they performed required 
parameter calculation uniquely, and they applied the technique of them to a Si0 2 system. An outline of 

the technique of Rappe and others is described below. 

[0223]Energy Ei by an individual electric charge of the atom i (Q) is considered to the 2nd order of 

Taylor's expansion, and is approximated with a following formula (1651). 

[0224] 

[Equation 84] 

Here, Xi 0 is electronegativity and it is expressed with a following formula (1652). 
[0225] 

Xj°=1/2 (IP+EA) — (1652) 

Jii is the coulomb repulsive force between the electrons on the orbit of the atom I, and is defined by 

the following formula (1661). 

[0226] 

Jii=IP-EA — (1661) 

As for IP, ionization potential and EA show an electroaffinity here. Total energy E (Q1, — , QN) of the 
molecule which comprises N atoms is the energy of an individual atom, and the sum of interaction 
electrostatic energy JyQi Qj between atoms. However, J- is a Coulomb interaction of the unit charge 

which exists at the center of atomic y. Therefore, a following formula (1662) is materialized. 

[0227] 

[Equation 85] 

The chemical potential of an atomic scale is given with a following formula (1663). 
[0228] 

[Equation 86] 

The conditions of the electric charge balance by this formula (1663), X1 = — If it substitutes for =XN, a 

following formula (1671) will be obtained. 

[0229] 

[Equation 87] 

That is, a following formula (1672) is obtained. 
[0230] 

[Equation 88] 

The conditions of all the electric charges shown in a following formula (1673) are taken into 

consideration. 

[0231] 

[Equation 89] 

Jjj is known and is a following formula (1674) about. 
[0232] 

[Equation 90] 

R is interatomic distance here at the time of an atomic size and I!=j at the time of I=j. Therefore, 
introduction of the procession C and the vector d will form the simultaneous equations showing in the 
following formula (1681) about Q. namely[Equation 91] 
It is CQ=-d if a definition is given. — (1682) 



It becomes. Although it was with C1 i=Qi and Cd =d in the paper, this invention persons used these, 
having corrected them as mentioned above. Required parameters, such as a size of atoms, such as Si- 
O-H, were computed uniquely, simultaneous equations were solved, and electric charge distribution was 
searched for. by this, a peak wave number is markedly alike and in agreement with survey. [ came ] 
[0233]An example of "a family register computed by the molecular dynamics simulator II", i.e., a vast 
quantity of data in which bond length (position vector r (t)) and a bond angle for every unit time are 
shown, is shown in drawing 24. This data is inputted into the membrane-characteristics calculation unit 
III, and basic membrane characteristics (optical, electric, crystallographic characterization) are 
computed here. 

[0234]About the process variation calculation unit I, an example of the two-dimensional general- 
purpose system is indicated to JP,H3-164904,A. 

[0235]As for a three-dimensional process simulator, this artificer succeeded in construction this time. 
Below, it has drawing 8 and this is explained without moving from its seat. Art and structure of a 
simulator are surveyed first. Although an application concerned has put a subject on an atom level, it 
has also proposed a new uniting method with other process variation calculation units I, is effective and 
more exact, and also shows that the new function which the former has can also be provided. A process 
variation calculation unit (10) of a partner who unites is the three-dimensional process simulator Ic, for 
example, and belongs to what is called a fluid model. In particular, this invention persons describe below 
the contents of the three-dimensional process simulator Ic who built newly. The main point is those 
with two and connection with an atom level simulator be possible for a "three-dimensional process 
simulator" whom this invention persons built first. Next, into a portion of this "three-dimensional 
process simulator", this invention person etc. are making full use of a chemicals model conquered newly 
and a mathematical expression. By these, union with an atom level was attained for the first time. A 
viscoelasticity coefficient, a stress constant, a diffusing constant, Cij, etc. are one of those are 
outputted from an atom level simulator, and these are transmitted to a fluid general-purpose three- 
dimensional process simulator. In drawing 8, it is (23), (33), (43), (53), (63), (73), and (83). For example, 
diffusing constants, such as an oxidizer, are transmitted to (23). A method of transmission procedure is 
mentioned later. Conversely, from a fluid general-purpose three-dimensional process simulator, it is 
transmitted to an atom level simulator and these viscoelasticity coefficients, a stress constant, a 
diffusing constant, Cij, etc. can be complemented mutually. That is, more precise prediction is realized 
by calculating by both three dimensions and exchanging a three-dimensional result mutually. 
[0236]By the general-purpose three-dimensional process simulator (10) side, calculation of process 
induction stress of an extensive field of several micron level according to input data by (2) and (3). More 
more exact functions can also be used by once transmitting a deed and this result to the molecular 
dynamics simulator II, or returning them to a general-purpose simulator (10) again via a file (22) and a 
file (32), and repeating this. Transmit a result of process induction stress calculation of a a certain 
grade extensive field which is a certain time, for example and for which it asked in a process variation 
calculation unit as a calculation procedure to the molecular dynamics simulator II, and once here, 
calculating a correction amount of a physical property value, anew, a motion of an atom is seen, a 
corrected physical property value is again sent into a general-purpose simulator, and results, such as 
process induction stress of an extensive field, are advanced by a general-purpose simulator, if 
sufficiently big stress is locally generated by repeating this, a motion of an atom can be seen by the 
molecular dynamics simulator II side, and propriety of plastic deformation can also be seen. 
[0237]Of course at the time of impurity redistribution calculation in the above-mentioned general- 
purpose process simulator portion, and stress calculation, a point defect action calculation part (9) 
described in drawing 8 is called. As a point defect, a vacant lattice in a Si substrate and Si between 
lattices are treated. Unless it takes in an action of these point defects correctly for diffusion of an 
impurity, redistribution of an impurity is not calculated correctly. Of course in an ion-implantation 
calculation part (4), a lot of vacant lattices by pouring damage and Si between lattices are generated. 
Although the redistribution of these is carried out to midium like a heat process, an impurity must have 
an interaction and they must make it reflect in calculation of impurity distribution then. Of course also 
by silicide, shape, and a stress calculation part (7), a vacant lattice in a Si substrate and an action of Si 
between lattices must be considered enough. 

[0238]by the way, drawing 8 — in addition to this, (8) is described for a few. At an LSI element trial 



production process, as for oxidation, diffusion, an ion implantation, etc., each required physical property 
value is already stored in this process variation calculation unit (10) based on many data. For example, 
PB1 of a formula (1 130) behind hung up about oxidation of Si and PB2 grade are it (1 131). D' etc. of a 
formula (1531) which hangs up a diffusing constant of an oxidizer in Si02 film behind are it. Generally, 
this appearance has data in general. 

[0239]However, data etc. are not necessarily assembled in LSI development. For example, trial 
production development of the memory cell using a ferroelectrics film is considered now. A case where 
PbTiO^ is formed is considered using CDV or a weld slag process. In the present general-purpose 

process simulator, yet, a diffusion coefficient of various impurities in inside of a PbTi0 3 film does not 

understand many so much, either, and the stress constant does not understand them so much, either. A 
segregation constant of an impurity in an interface is not known so much, either. In such a case, as it 
goes into a process of (8) first and is immediately shown in (84), it once goes to the molecular dynamics 
simulator II, and potential of a required element, etc. are computed, and a molecular dynamics portion 
moves there based on it, and a physical property value is computed. Although an already described 
portion is resembled very well after this, these physical property values are stored in a file from (83). 
And a calculation part is advanced like a process variation calculation unit. 

[0240]Here, explanation adapted to concrete processing is given with reference to drawing 9. Here, an 
oxidation process is taken up. That is, a wafer is thrown in, this is moved to a substrate treating 
process, and it goes into after an appropriate time at an oxidation process. 

[0241] Quantity of high purity oxygen gas is set to Pox, and distribution of this consists of this oxidation 
process as sigmatox, for example. ** and a temperature-up sequence are set to Temp, and constitute 
sigmatemp and a BARUPU sequence for distribution of this like tox and sigmatox. When having opened 
pulp, for example, a question is especially in a BARUPU sequence very here, but it is what was 
described in working example for details. 

[0242]For example, in drawing 8, oxidation, shape, and a stress calculation part of (2) are accessed, and 
the partial pressure Pox, the temperature Temp, etc. are inputted here. And stress calculation is also 
advanced to oxidation growth calculation and coincidence. At this time, distribution / change calculation 
is also performed using (21). Here, sigmapox and sigmate ** input and a variations from a center value 
is calculated efficiently. 

[0243]As for shape of an amount of growth of an oxide film, and an interface of an oxide film and Si, and 
also impurity quantity, these calculation results are calculated, for example. According to this invention 
person, stress is also called for. It is similarly from distribution and change. Shape of an amount of 
growth of an oxide film and an interface of an oxide film and Si and also impurity quantity are calculated 
also about Pox=Pox+sigma pox and portions, such as Temp=Temp lOsigmatemp. And these are once 
stored in a memory or a file (22), and are sent to a molecular dynamics portion. In a molecular dynamics 
portion, these are inputted as the signal a and position vector r (t) is computed. 

[0244]Although drawing 10 shows the same distance figure as drawing 9, it is related with formation of a 
ferroelectrics film here. 

[0245]By the way, explanation about a diagnostic portion in this invention is also given here. A part of 
binary tree by this example is shown in drawing 1 1. Here, a binary tree which paid its attention to an 
oxidation process is shown. As described also into a figure, Si02 thickness sent from the in-situ MO 2 
evening is first compared with Si02 thickness predicted from calculation here. And it is taking up and 
arguing about a time of thickness being settled in tolerance level. If it does so, 1R SUBEKUTORUKAPU 
sent from an inlitu monitor is compared with 1 R curve predicted from calculation as a following 
procedure. And as shown in a figure, first, a spectrum is compared and, moreover, a case where a 
spectrum shifts to the number side of Takanami is taken up by a case where it has shifted for a while. If 
it does so, it is possible as a following procedure that remaining stress has induced greatly in Si02. 
Then, the temperature can be slowly lowered in this case and remaining stress can be made to reduce. 
Then, an annealing portion indicated to drawing 9 is corrected, and it is inputted into (2) of drawing 8. 
On the other hand, when comparing a spectrum and going into tolerance level, it continues with the 
original process sequence. When thickness is not settled in tolerance level, poor processing is 
performed as correction being impossible. 

[0246]Data flow thinks in general that it understands now about drawing 8. I would like to make 



reference about an output of drawing 8 to a slight degree. That is, an executed result is once stored in 
a file (3). Here, a spectrum corresponding to a process searched for at each process, and optical 
[ other ] and electrical property results are stored. Although these give details later, they can be 
contrasted with quantity which is monitoring as data measuring. 

[0247]On the other hand, one person of this invention persons has proposed before a process simulator 
system which used variations calculation of a process variation. That is, they are JP.63-1 74331 ,A and 
JP,H3-164904,A. In these proposals, each process, such as oxidation/diffusion, an ion implantation, 
chemical deposition (CVD), etching, and lithography, can be dealt with by the handling like a fluid. For 
example, it is asking for redistribution of an impurity in a diffusion process by a fluid model, taking in an 
electric effect distinguishing an impurity total amount and electrical activation impurity quantity. 
[0248]An outline of a simulation system shown in JP,63-1 74331 ,A and JP,H3-164904,A is explained 
with reference to drawing 25 - drawing 49. In these proposals, it is treating as follows about a process 
variation. When a process of oxidation, an ion implantation, and lithography is taken for an example, for 
example, an average temperature (T) and with the permission Bala (deltaT), The permission Bala, a 
mean pressure (P) and with the permission Bala (deltaP) (deltat) and also an average mask size (L) and 
with the permission Bala (deltaL), an average dose (D), with the permission Bala (deltaD), etc. are input 
values. [ the average number of hours (t), ] 

[0249]As shown in drawing 25, a portion which calculates distributed change (deviation) to each 
subroutine, such as an ion implantation, oxidation/diffusion, and lithography, is attached to a process 
simulator. As a computational procedure, a method of calculating the variations and going is adopted, 
advancing calculation about these center values. 

[0250]This variations calculation is explained a little. If a temperature variations in a diffusion process is 
explained first, in a diffusion process, epsilon and the number of electrical activation impurities are set 
to N for an internal field according a diffusion coefficient to D and impurity distribution, an equation 
based on the second law of Fick is solved, and the impurity atom concentration profile C is searched 
for. A procedure of variations calculation is shown using drawing 26. Here, let left-hand side in a figure 
be a normal process (Normal process). 

[0251]A following formula (0071) can express a nonlinear diffusion equation. 
[0252] 

[Equation 92] 

k and k+1 are time here, and impurity concentration Ck+1 at the time of the time k+1 is an unknown. B 

of the right-hand side originates in the 2nd paragraph of a following formula (1072). 

[0253] 

**C/**t-XC k+r C k )/t — (1072) 

Here, since an implicit method is used, a coefficient matrix which appears in an upper type is altogether 
expressed using D, the mean temperature T, etc. It is diffusion-time tn+1 in the case of a normal 
process predetermined. Even eye watch is divided suitably and it goes in quest of a solution one by one 
to t1, t2, — , tn+1. Here, in each time, t1, t2, — , tn+1, these were linearized using the Newton method, 
and calculation is recommended. 

[0254]Next, a change process (Process deviated) is explained. Here, a case where change deltaT occurs 
to diffusion temperature is shown. When a variations of an upper type is taken and variation deltaCK of 
concentration in the time k is known, variations deltaCk+1 at a time of being the time k+1 is expressed 
by following formula (1 081 ). 
[0255] 

[Equation 93] 

That is, in the time k+1, the value D, psi, and N at the time of convergence of a normal process, C, etc. 
are saved, and F' is created using these, and the variation to the temperature T is prepared, and 
deltaCk+1 is calculated. Some matrix elements for creating F' here are shown in drawing 26. The hole 
model is fundamentally used for derivation of a diffusion coefficient. In this case, a Fermi level changes 
with temperature changes, the concentration of an electrification hole changes and a diffusion 
coefficient changes with these. Therefore, to the intrinsic carrier concentration ni which becomes the 
basis first, the variations comes to be shown in a following formula (1082). 
[0256] 



[Equation 94] 

deltanprijd.S/T+O.eOS/kT 2 ) deltaT — (1082) 

Change by the temperature of the internal field of an impurity can also be expressed using these. It took 
into consideration also to the variations of the increase phenomenon (Oxidation Enhanced Diffusion 
Effect) of the diffusion coefficient in an oxidizing atmosphere. In the oxidation process, since the 
linearity oxidation rate constant (Linear rate constant) and the parabola oxidation rate constant 
(paraboric rate constant) also changed, this was also taken into consideration. Previous F' is created 
using these. 

[0257]Here, a certain typical diffusion process is taken up and the situation of convergence in the one 
diffusion time is shown in drawing 27. The horizontal axis of drawing 27 shows number of occurrence, 
and the vertical axis shows the remainder. The converging state in a normal process is shown by 
drawing 27. In a normal process, it turns out that it is converging by 3 times of Newton loop. Variations 
procession F' previously shown using D at the time of this convergence, psi, N, etc. is created, and a 
temperature variations is calculated. As shown in drawing 27, it turns out that it converged by about 50 
times and has ended with 1/10 of the whole of computational complexity. 

[0258]Drawing 28, drawing 29, and drawing 30 show accuracy of a variations. Drawing 28 expresses 
impurity distribution at the time of 1000 ** over a two-dimensional lattice as a normal process. Drawing 
29 is a change process and is calculated based on drawing 28 by a case of delta T=1 **. On the other 
hand, drawing 30 hammers out 1001 ** as a normal process. 

[0259]value + drawing 29 of value = drawing 28 of these figures to drawing 30 — a value — it turns out 
that it is related. It turns out that it comes to be shown in drawing 31 when a value of drawing 30 is 
plotted anew, and it is in a tendency for impurity concentration to fall in the high concentration side if it 
shifts to an elevated temperature, and is in a tendency which increases in the low concentration side. 
[0260]When creating a semiconductor device, there are many processes and there are two or more 
parameters in each process. For example, an ion implantation to the things main at a CMOS process, 
and a (i)N-well field, (ii) It consists of patterning of Si 3 N 4 for an ion implantation to a P-well field, and a 

(iii)LOCOS process, a (iv)LOCOS process, a (v)n+ ion implantation, a (iv)p+ ion implantation, a (vii) 
annealing process, etc. When you consider fluctuation of a process, it is necessary to deal with change 
and frequency of many of these parameters (or frequency). 

[0261] Frequency of a variations and combination of execution are considered by JP, 63-1 74331 , A and 
JP,H3-164904,A. Now, since it is easy, three processes, alpha, beta, and gamma, are considered. In 
each, a center value of a process parameter C (Center), A value displaced to U (Upper) and negative in 
a just displaced value was named L (Lower), and from JP, 63-1 74331 , A and JP,H3-1 64904,A, those 
frequencies were defined, as shown in drawing 32. 

[0262]When the processes alpha and beta are considered, the combination can consider CU, UC, CC, 
LC, and combination of five pieces of CL, if calculation execution of UU, LL, etc. is disregarded. When 
passing through a process of gamma succeedingly furthermore, seven combination shown in a figure can 
be considered. However, he does not deal with these all but is trying to calculate from U in 
consideration of the symmetry of U and L in JP.63-1 74331 ,A and JP,H3-1 64904,A in L. 
[0263]Some trees of search of an inference engine which is the central part of this system are shown in 
drawing 33. Here, oxidation and a diffusion process are examined. "**" in a figure shows a predicted 
value by a statistical analysis system, and a "fruit" shows an actual measurement. 
[0264]An acceptable value used at this time is a value which passed through whole processes, such as 
oxidation, diffusion, and etching. At this time, epsilon value is obtained by performing a two-dimensional 
process statistical analysis system, when a "fruit" came out of "**" **epsilon, as shown in drawing 33, 
an oxidation process had abnormalities — purport description is carried out. In this case, temperature 
or time performs a load module of a statistical analysis system for which had shifted if needed. The 
above-mentioned assay is performed also about a lithography process, although explained focusing on 
oxidation and a diffusion process. Fluctuation by a lithography process is equivalent to a lateral slip of 
shape of a Si0 2 film, a junction type-like lateral slip, etc. 

[0265]Although Common Lisp is used for an inference engine, it is being suitable not only for numerical 
processing but processing of a text (List) and processing of shape, and because this is still very simpler 



for maintenance of a system. This system not only reads data of an actual measurement from an 
inference engine, but is performing a load module of a statistical analysis system. 
[0266]According to the simulation system shown in JP,63-1 74331 ,A and JP,H3-1 64904A A variation 
part can be processed now by computational complexity of the conventional abbreviation 1/100 by 
saving a coefficient of a procession at the time of convergence of normal repeated calculation in each 
process, and creating a primary variations procession using this, and taking into consideration change of 
each process parameter of said, and these crowdedness effects. By interlocking an inference engine 
and the above-mentioned system by Common Lisp, for example showed being made to a diagnostic pan 
of a process margin and also abnormalities in a process to optimization. 

[0267]ln JP,63-174331,A and JP,H3-164904,A, a basic diffusion equation is stood, as shown in a 

following formula (11 11) to each impurity. 

[0268] 

[Equation 95] 

Here, t shows time, x and z are coordinates of the horizontal direction of two-dimensional space, and 
the perpendicular (facing down) direction, and N is [ C is impurity concentration and ] electrical 
activation impurity quantity. D is a diffusion coefficient including concentration dependence and the 
enhanced diffusion effect, and a following formula (1112) can express it. 
[0269] 

[Equation 96] 

Di is a basic diffusion coefficient independent of concentration among a formula, and the temperature 

dependence is expressed by the following formula (1 1 13). 

[0270] 

D.=BBand (-BAAT) — (1113) 

The constant according [ BB and BA ] to the kind of impurity and k are Boltzmann constants. FE is 

given with a following formula (1 121) by making ni into intrinsic carrier concentration. 

[0271] 

[Equation 97] 

n is 1/2 C(C D -C A ) + (C D -C A ) 2 +4n 2 ]. Furthermore, it is ni, [Equation 98] 

ni =3.87x10 16 xT 1 - 5 x exp (-0.605AT) 

— (1122) 
It carried out. 

[0272]An oxide film growth rate was expressed as shown in a following formula (1 123) and (1 124) using 
the parabolic nature oxidation rate constant R1 and the linearity oxidation rate constant R2, 
respectively. 
[0273] 

R1=PB3and exp (PB1AT) — (1123) 
R2=PB4and ovn (PB2AT) — (1 124) 

When the above-mentioned parabolic nature oxidation rate constant is anew described as B and B/A, 
there is a relation which sets Zo and oxidation time to t and shows oxide film thicknesses in a following 
formula (1125) by making tau into a constant in respect of the horizontal coordinates x. 
[0274] 

These are prepared and it asks for a variations to temperature deltaT of each above-mentioned 
coefficient. 

[0275]A variations of n comes to be shown in a following formula (1131) as deltani =ni and 

(1.5/T+0.605AT2) deltaT. 

[0276] 

[Equation 99] 

A following formula (1 132) is materialized. 
[0277] 

[Equation 100] 

The variations of a growth coefficient is expressed by a following formula (1 133) and (1 134). 
[0278] 



[Equation 101] 

When it rewrites to B=R1, A=R1/R2, it comes to be shown in a following formula (1 141) and (1 142). 
[0279] 

[Equation 102] 

And the variations of t+tau becomes like a following formula (1 1 43). 
[0280] 

[Equation 103] 

The simulation system explained above was epoch-making as that time. It is impossible however, to 
already use for the process control and monitor art in the latest clustering tool. 
[0281]Then, this invention persons repeat examination further and came to complete this invention. 
Here, optimal design of the ferroelectrics film of the electric field effect type MIS element which used 
the perovskite single crystal was performed, the element was actually made as an experiment and 
evaluated according to this, and invention completion was checked. 

[0282]Namely, being provided by Claim 11 is in charge of forming an electric field effect type MIS 
element on a semiconductor single crystal substrate, In a system which uses a ferroelectrics film of this 
MIS element as a single crystal, and moreover contains this semiconductor substrate and this single 
crystal ferroelectrics film, It is a manufacturing method of an electric field effect type MIS 
semiconductor device considering it as design arrangement of a single crystal ferroelectrics film and a 
substrate which introduced a ******** valuation function into abinitio molecular dynamics theory 
expressing the object characteristic, and followed the optimal solution of a system under desired main 
movement conditions by this. 

[0283]Being provided by ** and Claim 12 is in charge of forming an electric field effect type MIS 
element on a semiconductor single crystal substrate, In a system which uses a ferroelectrics film of this 
MIS element as a single crystal, and moreover contains this semiconductor substrate and this single 
crystal ferroelectrics film, abinitio molecular dynamics theory expressing the object characteristic of the 
system — as a ******** valuation function — a system — it is a manufacturing method of an electric 
field effect type MIS semiconductor device arranging said single crystal ferroelectrics film so that the 
whole free energy may be taken up and free energy may become the minimum under main movement 
conditions. 

[0284]Being provided by Claim 13 on Si semiconductor single crystal substrate, In forming an electric 
field effect type MIS element, a gate insulation oxide film of this MIS element is used as a single crystal, 
and abinitio molecular dynamics theory which expresses the object characteristic of the system in a 
system containing this semiconductor substrate and this single crystal ferroelectrics film — as a 
******** valuation function — a system — the whole free energy, [ take up and ] This free energy is a 
manufacturing method of an electric field effect type SiMIS semiconductor device arranging said single 
crystal ferroelectrics film under main movement conditions so that it may become the minimum. 
[0285]Being provided by Claim 14 on Si semiconductor single crystal substrate, In forming an electric 
field effect type MIS element, a gate insulation oxide film of this MIS element is used as a single crystal, 
and abinitio molecular dynamics theory which expresses the object characteristic of the system in a 
system containing this semiconductor substrate and this single crystal ferroelectrics film — as a 
******** valuation function — a system — the whole free energy, [ take up and ] While this free 
energy arranges said single crystal ferroelectrics film so that it may become the minimum under main 
movement conditions, it is a manufacturing method of an electric field effect type SiMIS semiconductor 
device stopping oxygen deficiency concentration of that single crystal ferroelectrics film to 0.01% or 
less. 

[0286]Being provided by Claim 15 on Si semiconductor single crystal substrate, In forming an electric 
field effect type MIS element, a gate insulation oxide film of this MIS element is used as a single crystal, 
and abinitio molecular dynamics theory which expresses the object characteristic of the system in a 
system containing this semiconductor substrate and this single crystal ferroelectrics film — as a 
******** valuation function — a system — the whole free energy, [ take up and ] While free energy of 
this system arranges said single crystal ferroelectrics film so that it may become the minimum under 
main movement conditions, It is a manufacturing method of an electric field effect type SiMIS 
semiconductor device, wherein unevenness of local energy dedicated within deviation score 3sigma and 
stops oxygen deficiency concentration of the single crystal ferroelectrics film to 0.01% or less 



especially. 

[0287]These optimization procedures and a trial production result are described below. 
[0288]Here, as a single crystal ferroelectrics film, PZT system perovskite was taken for an example, for 
example, and a Si surface was made into a field (100). Cu of FCC structure may be sufficient as this 
substrates face as a metal electrode. The most important problem is transposed to a problem to what 
kind of crystal physical relationship single crystal ferroelectrics film PZT system perovskite and Si (100) 
side should be arranged. Here, total energy of a system was taken up as a valuation function which is a 
main point of this invention. And a Si/PZT interface was taken up as this system. Perovskite as PZT 
under impression conditions and physical relationship of Si (100) found out that it was optimal to 
perform it as follows by this invention persons' creating an energy expression as a strict valuation 
function beforehand, and using this. That is, the median line of an angle a1 axis of C4v structure 
expression of perovskite and a biaxial to make was aligned with [110] shaft orientations on Si (100) side, 
and c axis was leaned, and coinciding the 1st Si and the 3rd Si position with Si on Si (100) side found 
out that it was an optimum value. At this time, it was also simultaneously shown in PZT that less than 
0.01% of oxygen deficiency is tolerance level. A key map which this invention created to this appearance 
proposes was shown in drawing 34. 

[0289]It is attached to drawing 34 and explains in detail. An upper half of drawing 34 a shows a top view 
of a simple lattice by the side of Si, and a lower half is sketch drawing of a PZT side. It is indicated that 
drawing 34 c is best to arrange in parallel a line which ties perovskite to the direction of [110] Si. 
Although explained for details later at this time, in order for a straight line connected to silicon to 
become parallel to [110] Si, C axis of perovskite needs to incline. By carrying out like this, according to 
this invention persons' calculation, a line which connects Si to perovskite was only arranged in parallel, 
and it also found out that total energy fell 10 more% compared with a case where C axis is not leaned. 
When use for an element was further actually considered by besides and it was in PZT even within 
0.01%, it found out that the original feature of single crystal perovskite was shown enough. 
[0290]A problem of surface energy pointed out here affects an element characteristic plentifully so that 
it may explain below. An azygos atom will be made easy to form and energy by mismatching generated 
in an interface which a ferroelectrics film and a substrate form will promote formation of interface state 
density by this. This interface state density leads to hindrance of a conduction phenomenon, and a fall 
of reliability. A defect of an unusual bond etc. was made into a film and a possibility of becoming level is 
conceived. Here, tolerance level of these defects is shown for the first time. 

[0291]Comparison with the MIS characteristic and a conventional example which were created by this 
invention was also performed. A CV method was used especially here as what can evaluate interface 
state density effectively. Dispersion in data was also correctly shown in the characteristic by this 
invention. When based on this invention, it turns out that a value is falling to 10 by about 1/at the 
interface state density compared with a conventional example. When mobility of an MIS element by this 
invention was also evaluated under the same electric field situation, about 12% of improvement was 
found. Reliability also improved. 

[0292]A procedure which came to invent the axial new relation which this invention persons do not have 
until now is explained in more detail. Drawing 35 is structural drawing of completely new strict abinitio / 
molecular dynamics simulator which this invention person etc. created. An upper half of drawing 36 is a 
portion of abinitio, and a lower half is a portion of a molecular dynamics simulator, using this simulator - 
- for the first time — artificers — for example — Structure of perovskite showed what kind of thing it 
was. Drawing 37 is the first perspective diagram of perovskite in which this invention persons were able 
to show a strict structure. A large ball in a figure shows oxygen and a small ball shows Pb. Since it is a 
field complicated to this appearance which field of perovskite and which field of Si carry out **** 
junction as the Z-axis shows C axis of perovskite here and it understands by a diagram, it turns out 
extremely well that a guess is difficult. 

[0293]Paying attention to an oxide film, contrast work of spectrum data was done as an accuracy check 
of the above-mentioned molecular dynamics simulator. That is, ROKKIN mode, stretching mode, etc. of 
O-Si-0 or Si-O-Si were sometimes read in a position of Si atom of ****, or an oxygen atom, etc., and 
undamped natural frequency was calculated from these figures. When this is compared with an actual 
measurement, it is 490 cm-1 and 1 1 1 1 cm-1 to 440 cm-1 and 1 100 cm-1, and is an appropriate range, 
respectively. When a pressure was impressed, it also turned out that it is well in agreement with an 



actual measurement that intensity in O-Si-0 mode increases. 

[0294]This invention persons calculated also about a case where very big compression stress is 
impressed to PZT, also recommending such inspecting work. As an interface, a way which took such a 
structure also understood that there was a profit of total energy. 

[0295]A field etc. which were used for calculation are discussed in detail a little. In calculation, to the 
Si-substrate side, tens atomic layers were taken to a depth direction, tens atomic layers were taken to 
a depth direction in the surface, and tens atomic layers were taken also to a longitudinal direction, and a 
calculation area was taken. Also in perovskite, ten atomic layers of numbers took height in every 
direction at a time, respectively. 

[0296]It can be gathered by this method that new physical relationship was able to be found out by a 
technique twisted until now. This invention persons developed a strict addition type of strict ionic bond 
potential to a field containing a ferroelectrics film and a Si single crystal. In the conventional technique, 
stringency was missing, especially energy of a system was not able to be searched for correctly. In this 
Description, this invention person takes the side of a conventional fault and an inaccurate place, one by 
one, clarifies and goes. 

[0297]According to this invention person, it can explain as follows. This invention persons built the new 
simulator which is not in the former, and calculated about perovskite on a single crystal of a 
ferroelectrics film, and a concrete target as an example. 

[0298] Explanation is added a little, carrying out comparison with the former to below also about a 
simulator which this invention person created. For example, an oxide film is described first. Three kinds 
of potentials, between Si-O, between O-O, and between Si-Si, must be expressed strictly. Potential of 
total between Si-O, between O-O, and between Si-Si will be very difficult in practice, and all three 
sorts of amounts of addition will be infinitely emitted by a paragraph of the distance r. Since the 
Coulomb potential lengthened the skirt to a distance, the calculation itself could not close it, but it was 
using a method of Ewald by the conventional technique in somewhat strict calculation. However, 
although mentioned later, stringency is lacked in calculation. 

[0299]In actually arranging, a technique which this invention persons used is shown below. Here, this 
invention persons show the first simulator to create below. 

[0300]This invention persons built an original new abinitio molecular dynamics system, though it inquired 
wholeheartedly and computational physics was followed. It is why a new trial became possible just 
because there was this new system. What this invention person etc. created below is explained in detail. 

[0301]Although Si0 2 was taken for an example here, about other things, it also checked that it could be 
used only by replacing a number of an element, etc. 

[0302]An entire configuration of the new simulator which this invention person etc. proposed is 
described in drawing 38. The new technique of having compounded a molecular dynamics portion (MD) 
and a density functional (DFT) was proposed. Proper use with a molecular dynamics method and a 
density functional method by this invention persons is shown below first, using drawing 39. 
[0303]As shown in drawing 36, a molecular dynamics (MD) portion occupied a lower half of a flow chart, 
and a density functional portion (DFT) occupies an upper half, a molecular dynamics method does not 
come out OK — the so-called a little big system in limited temperature is treated, and potential 
beforehand searched for from the ultimate cause is used for potential between atoms (ion). On the 
other hand, a density functional (DFT) was made to use, when searching for a ground state of a small 
system. And naturally a portion of DFT set temperature to OK. It asked for a potential portion from a 
ground state of an electron system. A wave function (KS equation) was specifically solved and a partial 
density functional was used together. 

[0304]Calculation procedure of a simulator by this invention person etc. is actually explained taking the 
case of Si and O. It asks for total energy of a ground state of an electron system without degeneracy in 
a DFT portion, and is **. [External Character 1] 
It is]. 

[0305]namely[Equation 104] 
[External Character 2] 

N interaction energy and the 3rd paragraph are kinetic energy in case there is no interaction between 



electrons supposing electronic one-body approximation, [Equation 105] 

It is written. The 4th paragraph is exchange correlation energy and the form is based on partial density 
approximation. **[External Character 3] 
[Equation 106] 
[External Character 4] 

Equation (Kohn-Sham equation)[Equation 107] 
[External Character 5] 

It must solve self-contradicting [ no ] (self-consistant). 
[0306] 

[External Character 6] 

It applies and integrates with the electron density n to RUGIepsilon xc (n),[Equation 108] 

It carried out. The form of epsilon^ (n) is proposed variously the formula besides Wigner etc. Further, 

xc 

[External Character 7] 

A ** Coulomb field is not used as it is, but it sets by the norm preservation pseudopotential which 
erased the singularity in the nuclear position, and is a frog. The number of the ingredient in the Fourier 
expansion was able to be pressed down by accustoming in this way. 

[0307]As total energy, this invention persons use the thing also containing the coulomb potential 
between atoms (ion) in order to expect ****. namely[Equation 109] 
[External Character 8] 

epsilon etc. . A formula (2071) to the 1st paragraph of the right-hand side is the right from the formula 
(2065) of DFT.[External Character 9] 

It comes out. The portion is Veff from the formula (2065) of one-body approximation like the case of 
DFT.[External Character 10] 

A new view was put in here. It considers that E is potential energy and is kinetic energy virtual 
otherwise. [Equation 110] 

It introduces. Therefore, a Lagrangian is L=K-E. — (2073) 

Mj in a formula (2073) are mu and mu although it is the mass of a true atom (ion). It is virtual mass and 

changed according to the purpose like the after-mentioned. They are regular rectangular cross 
conditions as a binding condition. [Equation 111] 

Since it ****, when making the Eulerian equation, undetermined-multipliers lambda jk is introduced, and 
it is L.[Equation 112] 

The variations of ************ j s taken. The result[Equation 113] 

The temperature T of the kinetic energy of a formula (2073) corresponds. If it says by a relation with 
sigma1/2MjRj 2 , it will be a temperature that it is not virtual and actual, because the degree of means 

has temperature each one by the energy equipartition law of classical statistical mechanics, mu and mu 
nunu can be controlled or reinforced by size. 

[0308] For example, if it is referred to as T->0 from an elevated temperature as mu Mj, psij will change 

with Rj stopped and a ground state of an electron system will be acquired. The left side of a formula 

(2076) is set to 0 at this time, and the right-hand sides are a formula (2071) and the cautions under it. 
[Equation 1 14] 

Since it becomes, {lambda ik } which makes the linear combination of psij suitably is diagonal-line-ized, 

and obtains a KS equation (2064). That is, it means obtaining the same result as having solved the KS 
equation by DFT by simulated annealing. However, although the state of the temperature T is made, 
since it was not based on the Monte Carlo method but movement of the virtual dynamical system is 
followed, it is called dynamical simulatedannealing (DSA). 

[0309]When integrating with a formula (2076), psij cannot necessarily be treated as a variable as it is, 

but the expansion coefficient (that is, Fourier coefficient) which developed it by the plane wave can be 
treated as a variable. In order not to increase the number of this expansion coefficient recklessly, it is 
required like the case of DFT to erase the singularity in a nuclear position at pseudopotential. When 
especially the periodicity of a crystal can be used, the linear combination of psl which diagonaHine-izes 



above {lambda ik } is described as phi nR . k is a wave number vector in a Brillouin zone, and n is a zone 
index. 

[0310]phi nk [Equation 115] 

******** j s ** |ambda nk is a characteristic value of a procession {lambda ik } and is an energy level. phi nk 

(n, t) is developed. [Equation 1 16] 
[External Character 11] 

Effective potential V ff which is not filled is also developed.[Equation 1 1 7] 

A formula (2075) and a formula (2076) are substituted for a formula (2074).[Equation 118] 
[External Character 12] 

It becomes common, so remove. [Equation 1 19] 

****** Since this formula is clearing up the integration of the vibrating part analytically, purely, rather 
than a numerical method, it can enlarge deltat and can decrease computational complexity. 
[0311]The Coulomb potential is as follows, this invention person found out being divided into the 4th 
paragraph, when it decomposed. Although there was only the 3rd paragraph in particular conventionally, 
it checked that there was the 4th paragraph newly. These methods are described one by one. As shown 
below, this invention persons began to solve these as a start of a strict formula first from a basic 
equation which considered a dielectric constant. It is a basic equation first. 
[0312] 

[Equation 120] 

a conductor — epsilon=infinity and in the case of the vacuum epsilon= 1, the Coulomb potential differs, 
L is taken by 1 ** of a unit crystal (cube), sigma is taken within a unit crystal, and it is electrification 
and the position of Zi, and the oxygen i i~th particle. This is because polarization arises in the inner 
surface of a ball by electrification in a ball. Although the layer of electric doublet is made in the inner 
surface of the ball which is not a conductor, the last paragraph of the above-mentioned formula serves 
to negate the effect exactly. Since the method of Ewald gives the thing of the left side and 
epsilon=infinity, it must add the last paragraph of the above-mentioned formula in quest of the value 
within a vacuum. Only a result is hung up. The conventional derivation often has an error also in this 
portion. 

[0313]Coulomb potential[Equation 121] 
[External Character 13] 
It set up as follows. 
[0314] 

[Equation 122] 

It is expressed, n , n , and n 7 are the vectors of Tooru of x, y, and the direction of z here, respectively, 

x y z 

and n x , n , and n z are the integers ranging from -infinity to +infinity, respectively (in the case of a bulk 

crystal). It shall mean that ' of sigma' removes j-i at the time of n=o. 
[0315]New F function is introduced here, [Equation 123] 
[External Character 14] 
It ********(ed). 

[0316]This invention persons transform this still as follows, [Equation 124] 

Index m x , m y , and m z are the integers ranging from -infinity to infinity, respectively. = At the time 

[External Character 15] 
[Equation 125] 
[External Character 16] 
[Equation 126] 

[0317]The Coulomb potential which this begins, [Equation 127] 
[External Character 1 7] 
[Equation 128] 

This Since it does not contain, it is not related to power. It can change still like the next. 
[0318] 



[Equation 129] 

summarizing the result of former — coulomb potential — after all[Equation 130] 
[External Character 18] 

It is /L x , 0, and 0+m y (0 or 1"/L yf 0)+m z (0, 0 or 1-/L Z ). The good point of the condition of the above- 
mentioned expression is that decrease by the factor of erfc in phi (1), and the paragraph of sigma 
declines quickly by the factor of exp by phi (2) to the paragraph of sigma of a basis declining only to the 
order of a reciprocal in the paragraph of sigma. Since it is effective for reverse at slowness and 
fastness with attenuation by phi (3), suitable kappa which can balance both is chosen. It is the result of 
calculating contribution of Coulomb force sequentially from the near place of distance, and is a result in 
case the circumference is moreover a conductor, according to the case where the circumference is a 
vacuum — already — the 1st paragraph is added. Especially this is the portion which was not 
considered conventionally. 

[031 9]A point-of-rotation group which this invention person etc. created is shown below using drawing 
40. Si atom is put on the center (0, 0, 0) of a ball of the radius r, and four oxygen atoms are put on a 
position of the peak of a surface-of-a-sphere top regular tetrahedron. A position of C atom can be 
specified by an Eulerian angle (theta, phi). What is necessary is to face a vector (o, ogamma) which goes to 
the North Pole, and just to carry out rotation of the circumference theta of the y-axis, and rotation 
which is subsequently the circumferences phi of the z-axis first, in order to express it with coordinates 
(x, y, z). 

[0320]The rotation procession, [Equation 131] 

When one of four oxygen atoms is put on the North Pole (o, *), and moreover when [ here / phi is 
unfixed and ] *, the position of other three oxygen atoms, being undecided — an angle — psi — 
respectively — ( — 109 — degree — 28 — ' — psi — ) — ( — 109 — degree — 28 — ' — psi — + — 
120 — degree — ) — (— 109 — degree — 28 — ' — psi — + — 240 — degree — ) — expressing — 
having . Therefore, those (x, y, z) coordinates are formulas (2061), and costheta^-1/3, sintheta=2 root 
(2) / 3, and the thing that made phi psi, psi+120 degrees, and psi+240 degrees, respectively are made to 
act on t (0, Or). [Equation 1 32] 

It becomes. What is necessary is just to make R ( P hi, theta) of a formula (2061) act as it is, in order to 

make the 1st oxygen atom that had pointed out the North Pole direction of (theta, phi). As a result, this 

invention persons asked as follows. Coordinates (x, y, z) of four oxygen atoms, [Equation 133] 

It becomes. Coordinates (x, y, z) of four Si atoms in a unit cell are expressed with the parameter a, b, 

and c. When it compares with drawing 40, 2a is a cycle common to a lengthwise direction and a 

transverse direction. 

[0321] 

(I — 0, 0, and Oil III(a, b, c) (a-b, a+b, 2c) IV V (-b, a, 3c) — that is (0, 0, 4c), The 1st atom is set at the 
starting point, and if a z-coordinate goes up constant value every c and carries out orthogonal 
projection to a field (x, y), a square is made by I, II, III, and IV. A z-coordinate 4c Increases by the 1st 
atom of the following unit cell, and the 5th atom is in atom [ 1st ] right above. 90 degrees of 
arrangement relative to Si atom rotate at a time to a circumference of the z-axis, and arrangement of 
four oxygen atoms belonging to each Si atom goes. Movement of the above Si atoms is added to it. 
What is necessary is just to use coordinates (x, y, z) of four oxygen atoms of the circumference of the 
1st Si atom as it is, and a formula (2061) coordinates (x, y, z) of an oxygen atom around the 2nd, 3rd, 
and 4th Si atom, R (0, 90 ") rotates to an oxygen atom of the circumference of the 1st atom, 
respectively. Rotation of the parallel translation R of (a, b, and c) (0,180 degrees) Rotation of the parallel 
translation R of (a-b, a+b, and 2c) (0,270 degrees) It will be obtained if parallel translation of (-b, a, and 
3c) is made to act. 

[0322]In addition, although a group of I 111, III2, II III1 and IV2, and III IV1 and V2 is the respectively same 
oxygen atom, these give the same expression of relations as a top (it is natural). The oxygen atoms 11 
and 114 have shifted to a transverse direction by the one cycle 2a. 

[0323]When it had investigated to this appearance, it investigated about whether a process is further 
actually considered and it approves to defect density of how much. That is, now, for example, 3200 Si, 
and 3200 Ti were prepared, and 100 oxygen deficiencies were made from zero piece to this, this case — 
the next has sometimes comprised a position of Si atom of ****, or an oxygen atom, etc. clearly. First, 



relation of dynamic behavior of an oxygen atom and Ti atoms and a PZT structure factor was 
investigated. When an oxygen deficiency exists in the single crystal PZT, as compared with a case of 
being defect-free, turbulence is clearly received to distant Ti and oxygen of the Kazuhara child point, 
and, moreover, mobility of a point defect is a very large thing. It turned out that a direction in PZT is 
size very much about such dynamic behavior compared with a case of single crystal Si also in influence 
and a range. There is an ionic bond ingredient in PZT considerably, and characteristically it is easy to 
polarize in it, therefore a local deficit grows into a trigger, and such an active motion is also considered 
to be interlocked with potential change. This invention persons were such procedure, and when were 
investigated in detail, and they were 0.01% or less, they found out that an advantage of a single crystal 
ferroelectrics film could be pulled out enough. 

[0324] Bel ow is used about atom level calculation of a ferroelectric material by this invention person. 
This invention persons advanced PTO, structure grasp with an atom level of PZO and PZT, and grasp 
calculation with structure and membrane characteristics. A result of calculation is reported. As an 
object structure of using by atom level calculation, (1) With what has oxygen and Ti in the high position 
of symmetry with a perfect crystal of PTO, respectively, and a perfect crystal of (2) PTO. Oxygen is 
what (drawing 37) put an oxygen deficiency into one portion which met at c axis at thing shifted and the 
0.3A(3) above (2), a thing (drawing 43) which put an oxygen deficiency into (4) above (2) at one portion 
which met at an a-axis, the thing which replaced T of (5) PTO at Zr, etc. In particular, in (5), what 
replaced Zr in layers artificially was performed. It turned out that work of oxygen and Pb differs by PTO, 
PZO, etc. for the first time. From these, a cause of "spontaneous polarization" of a ferroelectrics film 
has been held. Relation of an atomic arrangement of a ferroelectrics film, distortion, deficit structure, 
etc. and basic membrane characteristics has also been held. 

[0325]Calculation with an atom level of the most fundamental PbTi0 3 film is shown first. About a case 
where hauling distortion is received, a case (drawing 44, drawing 45) where Pb position is able to be 
shifted from an equilibrium point for a trial, a case where it has an oxygen deficiency, etc., when 
searching for membrane characteristics, it was investigated whether electric charge distribution which is 
to foundations most could be found well. Drawing 44 is a picture of a unit cell of PbTi0 3 , and is a case 

where Pb is able to be further shifted a little from a center position. It asks for an isoconcentration 
figure of electric charge distribution at this time. Although drawing 45 is a unit cell too, it is a case 
where a portion of Ti-oxygen bond is made distorted. This result has shown the following points, (i) The 
whole has dramatically many ionic bond ingredients. It turned out that most electric charges can draw 
almost near to oxygen with strong (ii) electronegativity, and then the surroundings of Ti have an electric 
charge, (iii) — although electric charge distribution is symmetrical and spontaneous polarization cannot 
happen the way things stand, considering oxygen "deficit" and position changing of Ti, it turns out easily 
that local electroneutrality is broken. Furthermore by hauling, it has also checked that a dielectric 
constant fell. 

[0326]The following things were understood in a calculation result when it furthermore saw at 25 ** 
noting that there was no thermal agitation. That is, stable structure of PbTi03 ideal film which consists 
of the point group C4v was searched for. As a result, even if it shifted oxygen very only, it turned out 
that it returns to the point group C4v of origin with stable energy. Is there any more stable equilibrium 
point as it usually says also at a place which was still larger and shifted 0.3A of oxygen positions? This 
was investigated. It is becoming still more stable toward a position which is distant from 0.3 A of 
starting points very only. In a position usually said, this neighborhood is said for this gap, therefore 
spontaneous polarization to arise. It turned out that what has a stabilizing point in those with at least 
two and a cubic position with simple one, and another are the positions a little shifted for the first time. 
[0327]The origin of polarization is locally replaced also with an atom position which breaks down 
"neutrality" of charge quantity being a stabilizing point in free energy. In the above calculation, it turned 
out that existence of oxygen is very large. So, when carrying out an oxygen deficiency by accident, it 
turns out that it can become a big trigger also in crystal and electronically. (1) — it seems that 
localization of the electric charge distribution is carried out to three kinds each of atoms, respectively 
[ namely, ] 

[0328](2) According to an order in the periodic table O, and Ti and Pb, , respectively (2 or 2 s2 per s, 
two p4), it is (2 or 2 s2 per s, two p6, 3 s2, three p6, 3d2), and (2 or 2 s2 per s, two p6, 3 s2, three p6, 



10 or 4 s2, four p6, 14 or 5 s2, five p6, 10 or 6 s2 and six p2). A thing especially with a large electric 
charge distribution area of Pb is known reflecting this. 

[0329](3) When it sees well, divide and, as for electric charge distribution of Ti, it is in sight that the 
next O becomes irregular immediately. 

[0330]In order to deter this, it turned out that there is also a hand which replaces halogens, such as F, 
CI, Br, and I, well. 

[0331]example 2 In upper example, although described about Si-PZT, a spinel MgAI 2 0 4 film, cerium 

oxide Ce0 2 , and strontium titanate SrTi0 3 can also be used as a single crystal ferroelectrics film. 

Aluminum oxide aluminum 2 0 3 , zirconium oxide Zr0 2 , and yttrium oxide Y 2 0 3 , yttrium stabilization 

zirconium YSZ, Pr0 2 , and calcium fluoride CaF2 etc. can also use the cascade screen. It is possible to 

apply modification which sandwiches the silicon PZT also in structures other than this example to an 
interface with this single crystal ferroelectrics film, a ground silicon wafer, or a foundation electrode. 
[0332]If it is important to make free energy under voltage small and its attention is paid only to 
distorted energy, it is not necessary to be necessarily the minimum then. It is also based on thought 
that we showed perovskite. 

[0333]A actual system of measurement by this invention is also explained. Drawing 47 is a circuit of saw 
YATAUWA. A circuit of this method is explained. As for a high voltage transformer and Co, known fixed 
capacity and CRO of T are cathode ray oscilloscopes. If Co is chosen more greatly enough than 
capacity of a ferroelectric sample, most of secondary voltage of T will start a sample, but since it is 
suitably divided by the resistance Rd and is added to a horizontal axis of CRO, deflection of a horizontal 
axis of CRO is proportional to the electric field E in a sample. On the other hand, although deflection of 
a vertical axis of CRO expresses voltage between both terminals of Co, this is equal to DS/Co. That is, 
a figure of CRO expresses a D-E curve and an absolute value of D and E is also easily obtained from 
deflection of a vertical axis. The spontaneous displacement Ds is obtained as D when outside looks for 
a saturation portion to E= 0, and the spontaneous polarization Ps is decided by Ps=Ds. If paraelectrics 
without a loss are observed in this circuit, it will become a straight line which passes along the starting 
point, but a case where there is a loss, and a ferroelectric have hysteresis. A D-E curve which becomes 
this invention was expressed with a solid line to drawing 48, and a dashed line showed a result of the 
conventional PZT to it for comparison. What becomes this invention is understood that a dielectric 
constant is also dramatically large so that it may understand from now on. 

[0334]It is a place pointed out here at the impression time, and it tends to make the minimum energy by 
mismatching generated in an interface which a ferroelectrics film and a substrate form. And a thing on 
use was actually considered on this condition, and tolerance level of that quality was also shown. 
Reducing this surface energy leads to formation deterrence of level which becomes the hindrance of a 
conduction phenomenon, and improvement in reliability. According to this invention for this structure, it 
predicts by new abinitio / molecular dynamics simulator. Since any experiential model is not used for 
abinitio / molecular dynamics simulator which carries a desired valuation function, all places pointed out 
to this system can predict a natural phenomenon thoroughly. Here, also about a design, although an 
example which paid its attention to total energy at the time of impression of a system is shown, if it is a 
problem of a request of degradation deterrence, a design suitable for it will be attained further that a 
dielectric constant should be made a desired value, for example. 

[0335]How to calculate the physical properties of the two-dimensional above was shown. An example of 
a calculation method of three-dimensional physical properties was already explained with drawing 8. 
[0336]It is likely to be closed how **, and actual contrast and diagnosis are performed. That is, a 
simulation result is transmitted to a scheduler (13). At this time, from another side and a actual 
experiment system, monitor data is sent every moment and it comes. And data transfer from a file (4, 3) 
is synchronized with a thing with a scheduler as if these are stored in a file (4). And it is transmitted to 
contrast and assay / diagnostic block in the contrast diagnostic assay part IV. 
[0337]In a clustering tool used for this invention, drawing 49 is a key map showing how an in-situ 
monitor measuring device is added. As shown in drawing 49, for example, it is input gas composition 
analysis, a XPS signal, a FT-IR signal, etc., these are the average value in predetermined time space, 
and the amounts of monitors are a part of mere information. This invention persons provide a method of 



making a ************** measurand and sufficient information acquired from another side calculation 
contrasting. These are displayed as a picture [ graphics / on a monitor ] like drawing 50. 
[0338]Drawing 51 shows a temperature sequence as an example of input data to an abinitio molecular 
dynamics simulation system by this invention. Drawing 52 shows a concept of atom level execution 
prediction by an abinitio molecular dynamics simulation system by this invention. Thus, a reaction on a 
substrate, dissociation of a molecule, etc. are describing, for example with an atom level. 
[0339]Drawing 53 is a visualization figure of atom level execution prediction obtained by an abinitio 
molecular dynamics simulation system by this invention. After amorphous Si deposits this, it sees the 
state where it annealed at an about 600 ** elevated temperature, with an atom level. 
[0340]An abinitio molecular dynamics simulation system of this invention is applicable also not only to Si 
but an oxide film. That is, drawing 54 shows an example of execution about Si0 2 of an abinitio molecular 

dynamics simulation system by this invention. Here, a Si-O-Si angle in Si0 2 is called for every moment, 

and a spectrum of FT-IR is searched for only by calculation by carrying out the Fourier transform, 
asking for character frequency from now on, and superimposing this, using this prediction result — an 
actual measurement — what is necessary is just to carry out comparison assay with data every 
moment 

[0341]Drawing 55 shows an example of execution in an atom level using a simulation system of this 
invention at the time of using input data of drawing 51 shown previously. Drawing 55 shows performing 
as expected. 

[0342]From drawing 56, even if, as for drawing 69, annealing from membrane formation of a Si film at the 
time of using input data of drawing 51, distance between Si-Si which can be further set by solid phase 
growth, and a bond angle change time and temperature, it stops. 

[0343]In a system shown in drawing 1, basic membrane characteristics were computed in the 
membrane-characteristics calculation unit III. However, like a system shown in drawing 70, the special 
membrane-characteristics calculation unit III may be omitted, and the required characteristic may be 
computed inside process variation calculation unit I. They are deformable suitably according to a 
request of a design. 

[0344]As mentioned above, it became possible to reduce hardware and computational complexity 
substantially by this invention as compared with the former. For example, speaking of computational 
complexity, by the process simulator TOPAZ, as shown in drawing 98, a memory is about 24 M bytes 
and execution time becomes for several minutes. On the other hand, in process simulator SUPREM4 
[ famous ], since viscoelasticity etc. are taken into consideration, a memory amounts to 200 M bytes or 
more, and execution time will reach with clay also on the 1st. If it furthermore becomes molecular 
dynamics and ab-initio calculation, a memory has reached several times also several times and 
execution time. 

[0345]In such a situation, molecular dynamics and an ab-initio tool will be used for full, and prediction of 
atomic dynamic behavior and prediction of optical various characteristics will be started [ whether cost 
of a computer is large, and ]. This invention persons also find out the new use art of a simulator by this 
art while finding out a method of making full use of mathematical new art, and conquering these 
disadvantages. 

[0346]Conventionally, there are some which added some procedure above. It is a flow plan shown in 
drawing 99, and shows a main part of a flow of a density functional in the conventional molecular 
dynamics simulator. That is, in this technique, virtual mass is given and a Lagrangian of the whole 
system is used. This corresponds to a formula (1223) which this invention person etc. described 
previously. This is equivalent to the Car-Parrinello method already known. It is supposed that it is 
calculable of this method at a high speed apparently. However, a supercomputer of for example, a clay 
company level or calculation in one also still takes one whole day and night. Therefore, if 1-time 
calculation [ 1 time ] was done in order to investigate change of various parameters, it cannot become a 
diagnostic tool in real time at all. 

[0347]In order to actually work a process unit, appropriate preparation and expense start. According to 
this invention, without building a semiconductor manufacturing system and actually advancing a creation 
process of a semiconductor, only a process simulation or a thought experiment can be conducted and 
this is dramatically useful. For example, it becomes possible about looking for a new material easily and 



effectively moreover by conducting a thought experiment using various molecular structure models, and 
predicting the characteristic of a device using it. 

[0348]Drawing 131 is such a process simulator's lineblock diagram, and An input means of the internal 
memory means 3, such as the computer 1 and a hard disk drive, the mouse 5, the keyboard 6, etc., It 
consists of external storage of an output means of the monitor 7, the printer 9, etc., and CD-ROM1 1 
and floppy disk 13 grade, and its drive drive 15. 

[0349]A program which described procedure for performing a molecular dynamics simulation beforehand 
is created and stored in the internal memory means 3. This program is created by high-level languages, 
such as FORTRAN and the C language, combining assembler language as occasion demands, for 
example. Of course, such a program saved beforehand at external storage, such as CD-ROM and a 
floppy disk, may be suitably installed in the internal memory means 3. 

[0350]Although processing performed here is the same as what was already explained, Interface Division 
with a semiconductor manufacturing system which is actually working is omitted. That is, only 
processing corresponding to a loop including the molecular dynamics simulator II of drawing 1, the file 1, 
the process variation calculation unit I that calculates broad view character which is not taken into 
consideration to an atom level, and the file 2 is performed. 

[0351]Such a system can start and use a corresponding point from a program of a practical 
semiconductor manufacturing system of drawing 1. After realizing and putting a process simulator of 
only such calculation in practical use on a computer conversely, it is also possible to perform 
combination with a actual process unit. Namely, a semiconductor manufacturing system of drawing 1 
can be built from a system of process unit V, spectrum meter VI, and Drawing 131, An external output 
terminal of spectrum meter VI is connected to a communication port of the computer 1 using the 
suitable cable 19, and comparison processing of an actual measurement and a calculated value is 
performed. The result is sent to the control means 23 of a process unit similarly connected to a 
communication port of the computer 1 using the suitable cable 21, and is reflected in a process control 
in real time. 

[0352]About process induction stress, such as oxidation and silicide, a two-dimensional simulator has 
already reached even a utilization level mostly. However, this artificer discovered that a three dimension 
is indispensable for a design of a fine element, and completed this invention as a three-dimensional 
simulation so that it might explain below. 

[0353]In an old simulation, even if it uses oxidation and uses diffusion, a place of a join office belongs to 
a diffusion method rule of Fick. ** which that a substance does not move but situations differ greatly 
divides in a two-dimensional problem treating a x-y side when it is good noting that there is no 
movement of a substance in the direction of z which is the 3rd axis, but it is "stress." ****** and 
SHIRISAI Dison of this are the same as that of oxidation. 

[0354]Namely, stress calculation in two dimensions of TSUPEM4 said to progress most in the actual 
condition set stress to sigma, set displacement to v, viscosity was set to mu, and it set a BOASON 
ratio to nu, and was based on the following formulas. 
[0355] 

[Equation 134] 

sigmaxx+sigmayy=mu/(0.5~nu) - (**vxNo**x+**vy/**y) 
sigma xx-sigma yy = 2micro- (**vx no **x-**vy/**y) 
sigmaxy = 2micro- (**vyNo**x+**vx/**y) 

Understanding now is the handling of two-dimensional simple xy side, and it is that z ingredient has not 
appeared here. Although it is not visible explicitly here, Young's modulus E is treated as scalar quantity. 
[0356]This artificer performed creation of a strict three-dimensional rigid matrix, and preliminary 
calculations of a strict plane direction operator. When only the result was described, in Si system, it 
turned out that it is as follows. namely[Equation 135] 

**** of sufficiently thick two-dimensional plane strain (plane strain) is first considered to a depth 
direction using this formula. That is, about the case of (epsilonx!=0, epsilony!=0, epsilonz=0, 
epsilonyz=0, epsilonzx=0, epsilonxy!=0), when it substitutes for the above-mentioned formula, that of 
****** understands the paragraph of sigmaz. Namely, sigma z=Cxz-epsilon x+ Cxz-epsilon y+ Czi- 
epsilon xy. The value of sigmaz also traced that it was never slight quantity. Therefore, the first ** for 
that the paragraph of sigmaz which the three dimension is indispensable and should be taken into 



consideration also by two dimensions in the above-mentioned TSUPEM4 is missing to divide. Also in 
******, it became clear at least in two dimensions that z ingredient of the 3rd axis also comes out in 
the case of the stress of oxidation or silicide. 

[0357] Next, the above-mentioned formula was transformed and rigid-matrix Cij in each actual field was 

calculated. The elastic matrix of single crystal Si is described first. 

[0358] 

[Equation 136] 

Come out, and it is and is C1 1=1. 7 They are 2**1 5Gpa, C22=66**6GPa, and C44=-82**9GPa. 
Understanding now is never single scalar quantity. This is made to operate a field rotational operator, 
and it is as follows when it asks about the section (1 10) of the often seen [ first usually ] silicon crystal 
(refer to Drawing 132 (B)). 
[0359] 

[Equation 137] 

What is understood now is understood that contribution of z ingredient is very large. About the field top 

which becomes a problem by bipolar ** (111), it is as follows (refer to Drawing 132 (B)). 

[0360] 

[Equation 138] 

Although such converting operation was indispensable, in the old simulation, it was not taken into 
consideration but it became clear that the accuracy of the analysis of the fine structure is low. So, in 
this invention, a three dimension will be dealt with like the above. This is another aspect of this 
invention. 
[0361] 

[Embodiment of the Invention] Hereafter, various working example of this invention is described. 
[0362]working example 1 — this working example starts directly the prediction function of the optical 
spectrum which is one of the big features of this invention. It describes in detail once again here from 
the technique which this invention persons invented. 

[0363]First, a few is described from the conventional technical level about the art about a high 
reliability— ized oxide film design, its evaluation pertinent art, etc. As a publicly known publication, there 
is "ThePhysics and Technology of Amorphous Si0 2 " Plenum Publishing Corporation and ISBN 0-306- 

42929-2, for example. In this publication, the result of IR and the result of Raman are indicated as an 
optical measurement result of various oxide films. This publication is called thing outstanding especially 
as an argument of an atom level. 

[0364]It is argued by p. 63 of this publication by pole MAKUMIRIAN (Paul McMillan) how the Raman 
signal of an amorphous oxide film changes with processes. Drawing 71 and drawing 72 are the results of 
being measured by pole MAKUMIRIAN. For example, the curve a of drawing 71 shows the Raman 
measurement result of the usual oxide film, and the curve b shows the Raman measurement result of 
the densified oxide film. 

[0365]In the former, these both spectra should be compared, for example, moisture should **** by 
densifying, and the spectrum should change with these. Drawing 72 sees a spectrum before and after 
carrying out pressure impression for single crystal crystal. It is supposed that the place of 610 and 500 
cm-1 has change especially. However, it glances first, and the variation of a spectrum is very small and 
art is required very much about where its attention should be paid and whether it should be regarded as 
a capable difference so that it may understand now. A still more complicated experience is required for 
identifying the variation of a substance from the amount of spectral changes. 

[0366]Since it is in such a situation, in development of a new material which will be further needed from 
now on, there is a thing without the identification table of a spectrum plentifully, and the interpretation 
of a spectrum becomes a serious obstacle. 

[0367]From p. 71 of the above-mentioned publication, Hitoshi Kobayashi is seeing change of the Raman 
spectrum of the oxide fiber under high stress. A part of the result is shown in drawing 73. Here, the 
case where the case that stress is not impressed, and 3.5% of hauling strain are impressed is shown to 
the Si0 2 fiber. When 3.5% of hauling distortion is added so that drawing 73 may show, the place near the 

450 cm -1 has change. However, the variation of a spectrum is very small and art is required very much 
about where its attention should be paid too and whether it should be regarded as a capable difference 



so that it may understand also here. For identifying the variation of a substance from the variation of 
the spectrum, experience is required further. 

[0368]In the above-mentioned publication, Kobayashi has tried the interpretation of the spectrum using 
drawing 74. If a hauling strain is received so that drawing 74 may show, it is supposed that the 
tetrahedron portion of 109 or 5 times of Si0 2 structure will be distorted. However, he has only carried 

out only 1 set of tetrahedron structures, and is not arguing about what kind of relation it is as the next 
tetrahedron portion at all. Actually, since it is a certain limited temperature, it should certainly be 
accompanied by distribution of the angle. 

[0369]If quantitive distribution is not grasped as this invention person described on these Descriptions 
repeatedly, it will become about far from the analysis of a spectrum. The argument same as drawing 33 
is developed by Kobayashi. Drawing 33 (a) is Si04 as a function of an O-Si-O bond angle. It is a 
characteristic figure showing the wavelength of normal mode vibration of a tetrahedron. Drawing 33 (b) 
and (c) is a characteristic figure by the decoupling in the V4 mode produced by stress showing 
reduction in strong. 

[0370]Drawing 76 is a defect model figure in the oxide film by W.B. Fowler indicated to p. 107 of the 
above-mentioned publication, (a) of drawing 76 shows some coupled models of completeness alpha type 
Quartz, and are long combination and S of the inside L of a figure is short combination, (b) shows the 
model based on E1\ (c) shows an E4'model, and (d) shows an E2'model, respectively, (b) The unpaired 
electron in - (d) is expressed with e-. 

[0371]drawing 76 also shows — as — the above-mentioned publication — the tetrahedron structure of 
109.5 degrees — two pieces — or one piece is prepared and it is arguing about deficit structure. 
However, as mentioned above, it should certainly be actually accompanied by distribution of the angle at 
a certain limited temperature. Therefore, if calculation of a certain kind of many is not performed, it is 
far from the analysis of a spectrum about. 

[0372]Drawing 77 is an example of execution by the new molecular dynamics simulator which this 
invention person etc. created. Si used 32 pieces and oxygen used 64 pieces. There were 520 reciprocal- 
lattice points. A vertical axis shows the rate of change of volume, and a horizontal axis shows time, 
respectively. Becoming the external pressure and a balance mostly after [ of drawing 77 ] 1000 ft was 
confirmed. For example, the O-Si-0 angle, the Si-O-Si angle, and O-Si distance at this time are shown 
in drawing 24. 

[0373]The huge data equivalent to the "family register" of each atom is required so that drawing 24 
may show. This invention persons are preparing some devices also for the calculation technique. For 
example, what is necessary is to sometimes extract the mode of vibration with which the selection 
principle of infrared absorption is filled from the time series data of the vibrational motion of the atom of 
****, and just to ask for the frequency in infrared absorption. Anyway, such a trial does not still have an 
announcement and this can be said to be new "arms" of process development. 

[0374]Using a molecular dynamics system, enough, this invention persons set time as 0, and make even 
time t sampling time now after that in the place which reached the balance. Time correlation was 
searched for as follows, that is, it comes out. M (t) is polarizability and is expressed with a following 
formula here. Integral value of the above-mentioned formula - It is the value with which it integrated 
from infinity to infinity. 
[0375] 

M (t) =sigmae-r (t) — (1362) 

<M(t) -M (0)> expresses the autocorrelation function of M. 

[0376]ei is a charge number in the point here, and i is a number of particles. Although r is a vector and 
the measurement starting point does not have restrictions in particular, it is good to take on the square 
of calculation cells. 

[0377]What is necessary is just to ask for s (t1), s (t2), s (t3), s (t4), s (t5), s (t6), s (t7), and s (t8) in 
each time, in order to search for time correlation. 

[0378]for example, it crosses to the whole sampling time — sometimes — the calculated value of the 
coordinates position of ****, and its kinetic energy — quantity of motion etc. exist over all the particles 
further. And distance, those angular-position relations, etc. between still more nearly mutual particles 
are calculated. 



[0379]This invention persons calculated also about some crystalline oxide films and an amorphous oxide 
film. A part of the result is shown in drawing 78. That is, drawing 78 shows the calculation result of the 
Si-0 distance of single crystal cristobalite as one example. Drawing 79 shows the calculation result of a 
Si-O-Si angle in case there is a deficit. Drawing 79 shows around a deficit that movement is confused 
greatly. It also performed changing various concentration of these deficits and asking for the same 
action. 

[0380]Drawing 80 is an example of execution of the new molecular dynamics which this invention person 
etc. created, and simulates a motion of an atom. A large ball shows oxygen among a figure and the small 
ball shows Si. (a) of drawing 80 shows the atom in XYZ space, and (b) shows the projection to the XY 
side. Although it was not able to be shown to that locus here, in this example of a single crystal, the 
atom also mainly understood that the atom was exercising preferentially in the direction of the median 
line of XY axis. M(t) =sigmae-r (t) which this showed in the above-mentioned formula — (1371) 

A big bias will be shown in ****** and a spectrum peculiar to a single crystal will be shown as a result. 
[0381]This invention persons changed measurement temperature and performed the simulation also with 
the situation of the action of an individual atom, and the natural frequency based on this, and the case 
where distortion is added further. Drawing 81 and drawing 82 show those examples. That is, drawing 81 
is a figure showing a stress prediction result with same X-ray-spectrum prediction figure by the new 
molecular dynamics simulator by this invention and drawing 82. 

[0382]Two working example artificers enabled it to also perform prediction of the spectrum of an oxide 
film in case there is an oxygen deficiency, and enabled contrast with a actual oxide film. The example 
made computable to the optical property of an oxygen-deficiency oxide film as working example 2 is 
shown. 

[0383]When an ab-initio molecular dynamics method is used, it does not become a problem, but in the 
case of the molecular dynamics method which uses TTAT potential and BKS potential, when treating an 
oxygen deficiency, it becomes with a problem how best sets up charge of each atom. This is because 
combination is extended, or Si-Si combination arises and it becomes impossible to assume a fixed 
electric charge, when distortion becomes large dramatically. For example, if distribution of an electric 
charge changes actually and polarization becomes large by the elongation of Si-0 combination, it will be 
expected that the electrostatic energy of stretching vibration becomes large and the peak pitch of IR 
becomes large. 

[0384]Then, the balanced electric charge approximation (Charge equiliburation approach) using the 
ionization potential and the electroaffinity in which this invention persons had been tried by Rappe- 
Goddard for alkali and hydrogen, It applied to Si0 2 for the first time, and electrostatic potential was 

calculated by having incorporated amendment of the electric charge, and the technique of calculating 
movement and IR of each atom was performed. 

[0385]They are treated with balanced electric charge approximation, assuming the electrostatic 
energies of a system to be the ionization energy of each atom, an electroaffinity of each atomic pair, 
and a function of each atomic charge, It is the technique of solving a number equivalent to an atomic 
number of simultaneous equations, and asking for each atomic charge from conditions that the 
atomization study potential of each atom is equal, in an equilibrium situation. 

[0386]IR line spectrum and the conventional example at the time of using this method for beta-quartz 
were shown in drawing 83. By using this method, peak frequency shifted to the low frequency side about 
ten percent. This value is in agreement in survey and better accuracy. When a-Si0 2 also used this 

method, structure and an IR spectrum are enabled much more to relate with high degree of accuracy, 
and contrast with survey could be performed effectively. 

[0387]The key map showed intelligibly the method which this invention person etc. invented to drawing 
84. The search technique was shown in drawing 85. 

[0388]This invention is not limited to above-mentioned working example. Especially other devices, such 
as a liquid phase device and a superconducting device, are effectively applicable to manufacturing by a 
clustering manufacturing installation. In addition, it changes variously in the range which does not 
deviate from the gist of this invention, and is feasible. 

[0389]working example 3 — here, the details about the direct prediction function of an optical spectrum 



are shown as the 3rd working example. It creates on the computer of a-Si0 2 (amorphous Si0 2 ) first, 

and how to predict the optical spectrum is shown. a~Si0 2 was created by quenching, after fusing single 

crystal Si0 2 on a computer. Drawing 86 shows the temperature setting in the case of creation of a- 

Si0 2 . Although it may carry out slowly, melting skyrocketed temperature here and melting was carried 

out rapidly. Then, the temperature was lowered with the comparable temperature falling speed. Drawing 
87 is change of the pressure in that case, and drawing 88 is a volume change of a calculation area. In 
this process, this invention persons found out that the treatment of the external pressure and volume 
was important, when an experiential molecular dynamics method was used, or when an ab-initio 
molecular dynamics method was used and temperature was greatly fluctuated by rising and falling 
temperature, for example, in this example, if the external pressure is treated as a function of 
temperature so that volume may be maintained at about 1 law and a phase transformation may happen, 
and temperature rises, it will have set up so that a pressure may become large. The situation is shown 
in drawing 87 and 55. 

[0390]However, when attention is not paid to volume at all but the external pressure is kept constant, 
large expansion of volume as shown in drawing 89 arises, and it completely becomes actual elevated- 
temperature melting Si0 2 with Si0 2 of low density remarkably in things. Also by a cooling process, it is, 

even if it goes into solid phase, and Si0 2 of colander low density is formed from ****. Of course, when 

potential energy was investigated, it was dramatically high, and it was in the unstable state. If it 
becomes like this, the contrast with survey will be difficult. Therefore, paying attention enough to the 
external pressure and creating a~Si0 2 found out the desirable thing here, and it performed it. 

[0391]Change of the potential energy at the time of melting was shown in drawing 90. Since attention is 
paid to the external pressure, potential is rising gently, though it vibrates. Potential energy decreases 
gently in a similar manner also in the case of quenching. The temporal change of two sides of the field 
vertical to the cell c axis of the calculation area at this time was described in drawing 91. Since the 
technique of Parinello-Rahman is used for a cell so that shape may be changed according to a phase 
transformation, it is changing with melting. 

[0392]Since quartz is a rhombus at first, the angle of two sides to make is 60 degrees. However, since a 
crystal structure will break and internal pressure will be uniformly applied each neighborhood as for it if 
it fuses, the angle of two sides to make is approaching 90 degrees. It seemed that it did not return to 
quartz even if quenched, but a rectangular parallelepiped was maintained at about 90 degrees. 
[0393]When this invention persons created a~Si0 2 , they found out that there was a point which it 

should be careful of one more point. It is related with the potential of Si and O. Also by ab-initio or the 
molecular dynamics method of experiential potential (ATTA potential or BKS potential), in treating a 
Si0 2 system, In order change of the potential energy by change of atomic arrangement is dramatically 

large, therefore to treat atomic movement, I hear that it is necessary to ask for atomic movement with 
high precision than the algorithm of the usual Velret, and it is. 

[0394]Drawing 92 shows change of the potential energy at the time of asking for atomic movement with 
the algorithm of the usual Verlet, change of kinetic energy, and the sum of both energies. Total calories 
should become fixed when there is no exchange of energy with the external world. However, in the 
algorithm of the usual Verlet, since integration of the equation of motion drawn from raglan JAN is not 
performed in sufficient accuracy, a conservation of energy is not performed in this way. In order to raise 
integration accuracy, the situation did not change small the time interval which is an integration interval 
mostly at all. When this investigated, it turned out that it is peculiar to the material resulting from the 
potential of oxygen being dramatically steep. 

[0395]then — improving the algorithm of Verlet, in order to treat oxygen correctly — the time — 
****** — it asked until the individual atom of the acceleration of movement was high order, and it 
integrated with this, and saw in quest of speed and a position. Then, as shown in drawing 93, kinetic 
energy is changed by potential energy and an opposite phase, and total calories came to be kept 
constant, the exact treatment of such energy should do — when there was nothing, according to the 
result which this invention persons performed, at an elevated temperature, the badness of calculation 



precision broke in the case where atomic movement is intense, liquid Si0 2 , etc., structure broke 
disorderly by the cause, and creation of a-Si0 2 was not completed. The IR spectrum has also 
completely differed from survey. 

[0396]As mentioned above, this invention persons made it possible to create a-Si0 2 on a computer for 

the first time by elaborating uniquely in creation of a~Si0 2 . Thus, created a~Si0 2 had a little long Si-0 

bond length as compared with quartz, and, moreover, Si-0 length had fallen within the fixed range. In 
melting Si0 2 , although Si-0 of bond length smaller than the Si-0 length in quartz existed mostly, when 

it cooled, this of Si-0 length was lost, and he became only a combination longer than quartz. When the 
coordination number was counted and the first contiguity was set to 2A, 4 coordination was most, and 3 
and 5 were mixed for a while and they did not produce Si other than it. 2 coordination of oxygen was 
almost the case. 

[0397]Next, temperature was kept at 900 ** and structural relaxation of the time for 2 nanoseconds 
and the Si0 2 was carried out until potential energy reached a balance, in order to apply hydrostatic 

pressure or optically uniaxial stress to a~Si0 2 and to carry out prediction computation of the structural 

change of a-Si0 2 according to various stress fields. The structure of done various a-Si0 2 was changed 

into numerical values, such as each, an atom position, speed, acceleration, its primary differential value, 
a cell fundamental vector, and cell fundamental vector displacement speed, and was saved, respectively. 

[0398]After carrying out so far and creating various a~Si0 2 structures, it asked for the IR spectrum 

from each. It asked for the IR spectrum from the dipole moment. The dipole moments M are following 
formula M=sigma qi and ri. — (1461) 

************** is made. Here, qi is an electric charge of the atom i and ri is a position. When 
searching for this position and a system is neutral (sigmaqi =0), the value of a dipole moment becomes 
the same for anywhere about the starting point rO. It, [Equation 139] 

rO =sigmaqi, M=sigma qi ^(ri -rO) sigmaqi and ri-sigma qi, and ri -rO sigma qi=sigma qi, and ri — (1462) 
It is because it becomes. Since a paragraph which rO requires remains when a system is non-neutrality, 
an absolute value of a dipole moment changes with how to take the starting point, but when the Fourier 
transform is carried out, it turns out that frequency only appears in an ingredient of zero and does not 
influence others. In asking for this dipole moment, this invention persons have noticed that one-point 
cautions are required for how to take the atom position r. Since a periodic boundary condition is used 
for it when a molecular dynamics method is used for it, it is losing a discontinuous change of a dipole 
moment resulting from an atom left in - x direction of a cell entering from + x direction. 
[0399]ln consideration of being only carrying out a thermal oscillation and a-Si0 2 not being spread so 

greatly in temperature which actually measures IR, r used in the case of dipole moment calculation, 
Even if a periodic boundary condition was not considered but it went away from a cell in - x direction, a 
moment was calculated by having treated, when the position besides a cell had an electric charge. That 
is, it decided to deal with both of two atom positions, an atom position in consideration of a periodic 
boundary condition, and an atom position which does not take a periodic boundary condition into 
consideration, in the case of individual atom movement calculation. Thus, an example of a temporal 
change of a dipole moment for which it asked was shown in drawing 94. Since it was three-dimensional 
calculation, a dipole moment for x, y, and z all directions was obtained. 

[0400]In order to ask for an IR spectrum, as the first step, the Fourier transform of this dipole moment 
is carried out, and it asks for a power spectrum which is IR line spectrum. The example was shown in 
drawing 95. A peak position of an IR spectrum can be predicted in this stage. That is, it is that a big 
peak appears in about 1000 cm -1 , and that a big peak appears in about 500 cm" 1 . This result was in 
agreement with IR survey of a-Si0 2 . 

[0401 ]In order to ask for an IR spectrum, this line spectrum I (omega) is multiplied by a function of 

omega and the temperature T, and it asks by the following formula. 

[0402] 

[Equation 140] 



= [ alpha (omega) ] (4pi 2/3hcn) and omega-(1-exp (~h omega/kT)) -I (omega) 
— (1463) 

By calculating IR of a-Si0 2 in various stress fields, a stress field, the structure of a~Si0 2 , and the 

details (a peak shift, peak show RUDA, etc.) of the IR spectrum were associated. Thereby, contrast with 
survey was attained. 

[0403]The molecular dynamics simulator which this invention persons developed newly is applied to 
drawing 96 here at alpha~SiO 2 , The calculated IR spectrum is shown, the molecular dynamics simulator 

which this invention persons developed newly to drawing 65 is applied to a-Si0 2 , and the IR spectrum 

as which it was predicted when the external pressure was added is shown. The computational procedure 
by this invention persons is shown in drawing 100 and Drawing 101. A in the formula of movement of 
each atom within the block B in drawing 100 and B are paragraphs shown in a following formula. 
[0404]By using a system concerning working example 4 this invention, phenomena, such as a substrate 
and plastic deformation of an interlayer film, can be foreknown by a size of induction stress in process, 
prediction of the ingredient, and also the grasp of detailed relation with the stress and material physical 
properties. A stress problem is mentioned as an example especially in this working example. It was 
examined whether the reduction would be aimed at, having predicted stress induction how and taking a 
process variation range into consideration. 

[0405]It is indispensable to a fine element process design of these days to take into consideration 
stress at the time of process advance. About stress, the following phenomena arise, for example. 
[0406](1) In an embedding isolation method, trench a Si substrate and embed Si0 2 here. In this case, by 

the difference of a coefficient of thermal expansion between a Si substrate and Si0 2 , distortion and 

stress may remain and a crystal defect may be induced depending on a heat history. Therefore, it is 
necessary to draw up element structure and a manufacturing process which can reduce remaining 
stress as much as possible. 

[0407](2) It is becoming frequent to use a Si substrate with low interstitial oxygen concentration from 
the purpose of creating a quality oxide film. In such a situation, Si-substrate intensity decreases and it 
is easy to induce plastic deformation. 

[0408](3) For the purpose of formation maintenance of a shallow joining layer, and maintenance of a well 
area, it has increased as for a burden like heat processes, such as annealing and an oxidation process, 
in a substrate further by adoption of RTP. 

[0409]In such a situation, with how generating of a defect is deterred. If a value with larger stress than 
default value and its change predicted value is observed in the middle of a process, while clarifying the 
cause immediately, monitoring, in a future process, it is important that it is shown how correction should 
be added. This is much more important in a clustering tool. 

[041 0]A kind of stress and identification are very important further again. It is also important to also 
pursue principal stress and the principal surface of the maximum stress and to pay one's attention to 
decomposition shearing stress, although it is important. For example, when Si (100) board is used, the 
sliding surface of the highest priority is a field as shown in Drawing 102 (111). Although the [100] 
directions are slip directions on this (111) field, it is thought that a time of shearing force working is the 
most dangerous in this direction. 

[041 1]An output of these identifications is also possible for a simulation system by this invention 
persons. Thus, oxidation of a trench portion is first taken up for cautions with a purification to plastic 
deformation here, and a design manual for obtaining the shape of the oxidizing form, and stress 
distribution and further predetermined shape is described. For example, although it is said that stress 
has contributed "sharpening" and a "more round" phenomenon of oxidation of a trench angle, the 
actual condition is that the structure is not known yet. 

[0412]While understanding these phenomena in detail, it is required first to advance an exact process 
design which suppresses the above "sharpen" and advances "it is round." If induction stress is not 
grasped exactly, oxide film shape cannot be grasped correctly, either, a problem of stress is only a 
problem of plastic deformation — they are also not **** but shape grasp, and an item further 
indispensable also to impurity distribution grasp under a stress field. Although prediction of oxide film 
thicknesses after a LOCOS process of a SOI substrate which will be developed with an element after 



1G of future is divided for the time of oxidation advance unlike the usual Si substrate and it is easy to 
be influenced by stress, the quantitative prediction is a more complicated problem. 
[0413]In such a background, work must be extended instead of the conventional oxidation diffusion 
model even to action grasp of an oxidizer under stress, and an understanding of an oxidation 
phenomenon under stress and impurity re-diffusion under stress. And it must be calculated at high 
speed for a short time. This invention persons created generating and prediction of redistribution of 
stress, and the still more nearly computable new high-speed general-purpose process simulator portion 
to a part for a point defect two-dimensional diffused part under stress, and added this to an atom level 
simulator. That is, the high-speed two-dimensional oxidation / diffusion / stress / modification process 
simulator part by original development were added to a parent of an atom level materials design system. 

[0414]Survey of in situ observation of stress used micro Raman and FT-IR. A calibration of these 
measuring instruments was performed beforehand. When an operation sheet is inputted, Si and the 
amount of many [ physical-properties ] about Si0 2 are beforehand performed immediately by an atom 

level materials design system, and are calculated. This is sent to high-speed two-dimensional 
oxidation / diffusion / stress / modification process simulator part. 

[0415]This invention persons first describe phenomenon analysis of "sharpening" and "more round" 
oxidation traced using this new high-speed two-dimensional oxidation / diffusion / stress / modification 
process simulator part, and show a basic examining result of an embedding isolation process, etc. as an 
example of application further. Cubical expansion will be caused, if Si oxidizes and it becomes Si0 2 . 

Although there is viscoelasticity in Si0 2 and a part of that stress is eased at this time, stress still 

remains to both Si and Si0 2 . 

[0416]Oxidation of a trench corner is taken for an example, such a situation is investigated, and what 
was illustrated is Drawing 103. First, at this time, although an oxidizer of the concentration C reaches 
Si0 2 and invades and goes into a Si0 2 film, also when [ at which diffusion D' of an oxidizer is influenced 

by the stress sigma in a Si0 2 film ] it is a Si/Si0 2 interface and an oxidizer reacts to Si further, it is 

influenced by remains accumulation stress in a substrate. 

[0417]As a concrete effect, when the compression stress sigma exists in Si0 2 , a network atomic 

interval will control diffusion of narrowing and an oxidizer, for example. When the compression stress 
sigma is accumulated on a field of a Si substrate, Si atom interval will suppress invasion of narrowing 
and an oxidizer, and will stop an oxidation rate as a result also here. 

[041 8]If an ingredient of a x direction of the principal stress S is set to sigmaxx and an ingredient of a y 
direction is set to sigmayy now, about stress on a two-dimensional side, a following formula (1531) and 
an expression of relations - (1534) will be materialized, mu is coefficient of viscosity, nu is a Poisson's 
ratio here, and u is displacement. 
[0419] 

[Equation 141] 

By the way, generally the temperature dependence of the coefficient of viscosity mu can be expressed 
by the above-mentioned formula (1534). Here, VISC.O and VISC.E show a proportionality constant and 
activation energy, respectively. Therefore, when the above-mentioned angle stress type is seasoned 
with the temperature dependence of mu, it turns out that the ingredient of each principal stress 
decreases with a rise in heat. Reaction constant kappas ' in case the oxidizer which reached even the 
Si/Si0 2 interface reacts in a Si surface can be similarly expressed in writing like the above-mentioned 

formula (56). It is a value in case, as for kappas, stress does not exist here. Here, sigman is sigman =~ 
(sigmaxxnx 2+sigmayyny 2+2sigmaxynx ny), and is sigmatau=- (sigmaxxny 2+sigmayynx 2-2sigmaxynx 
ny). VR is a parameter. This formula shows that an oxidation reaction constant decreases, while 
compression stress increases. The influence of the remains accumulation stress in a substrate in case 
an oxidizer reacts to Si by the above-mentioned Si/Si0 2 interface is this. 

[0420] Diffusion coefficient D' of an oxidizer in an oxide film can be expressed by the above-mentioned 
formula (1531). p is expressed by following formula here. That is, it is p=-1 / 2 (sigma xx+sigma yy). D is 
a diffusion coefficient in case stress does not exist. Implications of this formula can be interpreted as 



follows. That is, stress sigmaxx and sigmayy take hauling to positive, and have taken compression to 
negative. Therefore, compression stress exists in a film, and if a substance is in a dense state, it is 
shown that a diffusion coefficient becomes small. The above-mentioned formula expresses this. 
[0421 ]A computational procedure of only this portion which this invention persons found out is shown in 
Drawing 104. First, it asks for concentration of an oxidizer in each node in a Si/Si0 2 interface. After an 

appropriate time, a temporary oxidation rate in each 1st node is computed based on this concentration. 
Moreover, the amount of growth u of a Si0 2 film in each node is calculated. It asks for the stress S 

from the above-mentioned formula (57) using this u. using this stress — a formula (56) and a formula 
(1531) — an oxidation rate is further recalculated again using a formula (1534). And again, an amount of 
growth is calculated and stress is computed. And if enough completed by stress, advance of new 
oxidation time will be begun anew. That is, a Newton's method is used. 

[0422]As shown also in this Drawing 104, stress calculation at the time of oxidation requires very 
prolonged calculation time. By typical trench structure, 900 ** oxidation was tried by 1300 initial nodes. 
An oxidizing atmosphere was made into 100% of dry oxygen, and oxidation time was made into 300 
minutes. Computation time was investigated using EWS Sun4. If it is conventional technology, time for 
15900 minutes will be required in general. This reaches also in about ten days and it is very unsuitable 
for practical use as [ this ]. 

[0423]About calculation of a Newton's method portion, this invention persons introduce a ** part 
principle, save the last converging state, and advanced convergence calculation using primary 
differentiation of this. By this method, 900 ** - 1 100 ** were able to be taken for an example, and 
conventional method is short to 1/1000 or less, and it was able to carry out. Thereby, this simulator 
advanced rationalization of a material property value, and improvement in the speed of operation part, 
and reached till a place which can be used as a tool for real time. This invention persons investigated 
accuracy of a time dependency value of oxide film thicknesses, the shape reproducibility of a trench 
corner, etc. as basic accuracy of a simulator further. 

[0424]Drawings 105-77 summarize a result of having investigated temperature of growth thickness, and 
time dependency in a graph using the above-mentioned simulator. Although not displayed here, oxygen 
tension questioned this invention persons also about a decompression action at the time of 10% or 1 
more%. As a result, since it goes on while stress accompanying expansion of an oxide film eases enough 
when an oxidation rate is slow, stress is not accumulated in a substrate, therefore a uniform oxide film 
is formed also in a convex portion. Oxygen tension was gradually understood that a method of 
mentioning and going is effective, dropping oxygen tension, making an oxidation rate late in early stages 
of a result of this invention persons' prudent examination, especially oxidation, and making it ease 
enough. 

[0425]This invention persons investigated why "it is round" would occur, saying "sharpen" of a trench 
angle using this system. As already stated, it is necessary to take some elements into consideration in 
oxidation of a trench angle. Five backbone corresponding to the following is described. 
[0426](1) Inside of these [ Si-surface stress-dependency model this invention persons of the plane 
direction dependency model (5) oxidizer / Si reaction of a stress-dependency model (4) oxidation rate 
of oxide film visco-elastic model (2) silicon flexible model (3) oxidizing agent diffusion ] five models, But 
a dominant thing traced that it was stress of a Si substrate of the 5th model, i.e., a Si/Si0 2 interface, 

also unexpectedly. Here, dry oxidation in 900 ** will be taken up. If shape of an oxide film at this time 
and shape of a Si/Si0 2 interface are often seen, shape "is sharpened" still will not almost be accepted. 

When it became after 300-minute progress by 900 ** dry oxidation, in oxidizing temperature of about 
900 **, a phenomenon in which a trench angle sharpened was analyzed as follows. 

[0427]That is, if it is till about 100 minutes, stress will not be generated so much. However, compression 
stress concentrates on the inside of a Si0 2 film of a corner, and the Si tip part side by Si0 2 cubical 

expansion with oxidation. At about 900 **, the viscous flow effect is small and a shape relaxation effect 
of an oxide film is especially small. And oxidation reaction is deterred by compression stress in a Si 
surface, and it sharpens as a result, however — if it becomes 1 100 ** — the viscous flow effect — 
**** — a shape relaxation effect of an oxide film which it hears sets, and is — ** And there is no 
accumulation of compression stress in a Si surface, and it becomes round as a result. Drawing 108 



showed this situation. The predicted value of these stress was able to carry out the measurement 
check in less than **10% of accuracy as a result of the Raman measurement or IR measurement. 
[0428]Preparing the above knowledge, this invention persons applied the above-mentioned simulation 
system to a part of basic examination of STI (Shallow Trenchlsolarion) for the purpose of area reduction 
of an isolation region further. Basic examination of a more round process of an upper corner of a trench 
was advanced especially here. Oxidizing temperature is 1 100 **, oxygen tension is 3%, and it saw about 
the 100-minute progress back. As a result, the following things are understood. 

[0429](1) In the beginning, oxidation of poly Si and a crevice between retreat parts of a buffer oxide film 
oxidize. At this time, a corner sharpens and is feeling. 

[0430](2) Oxidation progresses succeedingly, a crevice oxidizes, if filled, supply of an oxidizer will 
decrease rapidly, oxidation from the side will become rate-limiting, and a corner will be roundish as a 
result. The neighborhood of a crevice of principal stress is large. 

[0431 ]On the other hand, stress distribution at the time of oxidation by 50% of 1000 ** oxygen tension 
was examined. Time is after 30-minute progress. As a result, compared with stress distribution of 
previous 1 100 ** oxidation, it turns out that stress is large a little. The above consideration shows that 
it is the first condition not to store up remaining stress in a Si surface, in order to deter **** of a 
trench corner. It can consider stopping an oxidation rate enough and making modification relaxation of 
an oxide film cause enough also at low temperature as a fundamental measure for this etc. 
[0432]This invention persons advanced union-ization with a plastic deformation simulator further. As a 
result, it turned out that not only oxidation but modification after TEOS film deposition, etc. are 
investigated. 

[0433]Based on the above result, especially this invention persons saw a value projected on a field 
(1 10) with a two-dimensional simulation system paying attention to decomposition shearing stress on a 
field which is a sliding surface (111). If heat unevenness occurred within a wafer 300 degrees as a heat 
history although what is necessary is just to have been a part for /, and even a part for /was reached 
600 degrees selectively, it also traced with stress that a rearrangement occurred depending on a part. 
[0434]The above reserve knowledge is prepared separately and this invention was applied to a 
manufacturing process of a bipolar CMOS in real time. That is, contrast, and the amount of dispersion 
and a permissible dose of serial prediction by a simulation and the minimum necessary in each process 
which are a factor of each process further, a spot measurand, were examined. An example carried out 
about oxidation time dispersion of sacrifice oxidation and permissible doses, such as shape, is shown 
below especially here. 

[0435]A basic process is as being shown in Drawing 109. Furthermore it continues and there is each 
process of sacrifice oxidation -> of a RIE etching -> oxidation -> polycrystal Si embedding ->CMP-> 
polycrystal Si portion for isolation. A factor of each process has the range of fluctuation. In this working 
example, it is a portion of sacrifice oxidation of a polycrystal Si portion, especially its attention was paid 
to oxidation quantity. As the range of fluctuation of oxidation quantity, a real-time calculation result is 
shown about two, a case of oxidation quantity zero, and the specified quantity. Change of oxidation 
quantity, oxide film shape in a LOCOS process of continuing the shape of an oxidizing form of that way, 
stress distribution, and also next, and the ingredient characteristic of that stress is shown. 
[0436]Drawing 110 a distribution map of 1000 ** and sigmaxx at the time of an end of a LOCOS 
process, and Drawing 111, Similarly a distribution map of 1000 ** and sigmayy at the time of an end of a 
LOCOS process, and Drawing 112, Similarly a distribution map of 1000 ** and sigmaxy at the time of an 
end of a LOCOS process and Drawing 113 a distribution map of sigmaxx when temperature is lowered 
to 25 ** after an end of a LOCOS process, and Drawing 1 14, Similarly a distribution map of sigmayy 
when temperature is lowered to 25 **, and Drawing 115 show a distribution map of sigmaxy when 
temperature is lowered to 25 ** after an end of a LOCOS process after an end of a LOCOS process, 
respectively. 

[0437]As shown in these figures, when there is no sacrifice oxidation over a polycrystal silicon film, the 
surface of a polycrystal silicon film and the surface of a nitride have gathered after CMP. When a 
LOCOS process is performed in this state, an oxide film which expanded by oxidation will push from the 
side of a nitride. On the other hand, since a nitride bends and goes well when sacrifice oxidation is 
performed also with some, stress eases and goes. The result was shown in Drawing 116. 
[0438]This invention persons use the above which this invention persons have advocated, and a 



simulation system which made a high-speed molecular dynamics simulator using experiential interatomic 
potential, and a general-purpose process simulator unite, An oxidation process, an interlayer insulation 
film formation process, and after passing through a wiring process further, it traced that stress which 
remained to a Si substrate and an oxide film of a MOS device had big influence on the electrical 
property of an oxide film especially. 

[0439]While this invention persons predicted a stress generation process and a generation part in real 
time one by one using the above-mentioned simulator, for example, what kind of holding property a 
tunnel oxide film of an EEPROM element portion shows eventually made them output. An element in a 
present progressive can predict so much what kind of the characteristic is given to allowable width of 
the holding property of an element obtained eventually by using this system. This is the effect of a 
system which becomes this invention. A process of continuing can be corrected or avoided so that a 
system of skillful Lycium chinense may be used and the optimal characteristic can be obtained during 
process progress. This working example is shown below. 

[0440]working example 5 — this working example shows an example which applied this invention to a 
stress problem in a silicon oxidizing film. A stress problem is important in order to control a crystal 
defect in inside of a silicon substrate, and generating of a rearrangement but not only in it but various 
materials, it is important also as a factor which affects membraneous quality. 

[0441]In a silicon oxidizing film, oxidation stress resulting from cubical expansion at the time of silicon 
oxidizing and heat stress produced by rising and falling temperature have arisen, and these have 
affected quality of an oxide film. That is, if it says with an atom level, stress may affect network 
structure of silicon oxide and may influence even the Tetrahedron structure which is unit structures of 
a Si0 2 system. Influence is also brought to electric insulation if it lengthens. Then, this invention 

persons advanced an oxidation process, controlling quality of an oxide film, preparing structure in an IR 
spectrum and an atom level in silicon oxide, and knowledge to an electrical property. 
[0442]First, this invention persons investigated structure and an electrical property in an IR spectrum 
and an atom level of silicon oxide as a preceding paragraph story. As shown in Drawing 117, a silicon 
(100) board was prepared, an isolation process was performed, and an element region separated by 
isolation region was created. Next, a silicon substrate was introduced into an oxidation system, and as 
shown in Drawing 1 18, a gate oxidizing film was formed. At this time, oxidizing temperature was variously 
changed from 800 ** to 1 100 **. 

[0443] Substrate temperature was made to lower after forming silicon oxide of predetermined thickness. 
At this time, a temperature fall of substrate temperature was performed in two steps, as shown in 
Drawing 1 19. The 1st step is the process of lowering substrate temperature to 900 **, and made 
substrate temperature quench. The 2nd step is a process made to lower from 900 ** to a room 
temperature, and consideration in particular was not carried out to a temperature falling speed. When 
oxidizing temperature was set as 900 ** or less, the temperature was lowered only in the 2nd step. 
[0444]It opted for setting out of this temperature fall condition from analysis of a viscous relaxation 
action of an oxide film separately. Generally, at an elevated temperature, an oxide film becomes soft and 
is said to become easy to ease stress by viscous flow. Then, this invention persons investigated a 
stress relaxation action by viscosity in detail. An example which investigated change of the 
characteristic of stress relaxation by a temperature falling speed was shown in Drawing 120 and 
Drawing 121. In this example, oxidizing temperature was 1000 **. 

[0445]When it holds at 1000 **, respectively and quenches by 600 ** / min after Drawing 120 
performed a typical oxidation process, it is the result of investigating the maximum of principal stress of 
a direction along the surface at the time of cooling slowly by 6 ** / min. Stress eased in about 1 
minute, by quenching, viscous relaxation was slight and, a case of 1000 ** maintenance, and in the case 
of annealing, a result shown in Drawing 120 showed only being superimposed on heat stress, and that it 
was. 

[0446]Drawing 121 is a figure which made a horizontal axis temperature and plotted a stress action in 
the same temperature fall process. Stress relaxation is remarkable, an action of a viscous body or a 
viscoelastic body was shown, and stress relaxation hardly happened below 950 **, but Drawing 121 
showed that an oxide film showed an action of an elastic body in a not less than 950 ** temperature 
range. 



[0447]If the two-step temperature fall performed by this example is used so that it may guess from 
such example computation, By passing a temperature range which viscous relaxation produces as much 
as possible for a short time, it becomes possible to suppress viscous relaxation and various heat stress 
produced from a difference of a coefficient of thermal expansion of a silicon substrate can be set to an 
oxide film by moreover changing various oxidizing temperature. 

[0448]Thus, when a tensile stress of the direction of a substrate face directly under silicon oxide was 
measured from a section using a micro Raman spectrum using a created sample, it was checked that 
various heat stress has induced in a silicon substrate. This invention persons compared a calculation 
result illustrated to stress induction measured value and Drawing 122 in these substrates, and guessed 
a stress value induced in an oxide film. In performing these stress calculation, this invention persons 
used material and a process prediction simulation system which a molecular level to a general-purpose 
level which this invention persons built included. As a result, a relation shown in Drawing 123 between 
oxidizing temperature and compression stress in an oxide film was obtained. Drawing 123 shows that the 
maximum of compression stress in an oxide film increases linearly with a rise of oxidizing temperature. 
[0449]An IR spectrum was measured to an oxide film which produced compression stress formed as 
mentioned above as next preparation. On the other hand, computer prediction of oxide film structure 
when this compression stress is added, and IR PEKUTORU was performed using a classic molecular 
dynamics method. Drawing 124 shows change of calculation IR PEKUTORU by existence of compression 
stress addition. 

[0450]Drawing 125 shows the result of measuring a wave number shift by stress of a Si-O-Si 
asymmetric stretch peak from measured value of an IR spectrum. When an actual measurement was 
compared with a calculated value, signs that an appearance wave number of an asymmetric stretch 
peak of Si-O-Si shifts to the low wave number side rather than a case where he has no stress addition 
are obtained from a calculation IR spectrum with compression stress, and it turned out that an 
experiment and calculation are well in agreement. Then, amorphous structure and IR full spectrum at 
the time of adding stress variously were calculated and prepared. 

[0451] Furthermore, relation of stress and an electrical property was investigated as next preparation. 
An oxide film which has produced the above-mentioned various compression stress was used as gate 
dielectric film, a gate electrode was deposited on it, and Qbd was measured. Deposition of GE 1 ****** 
is 900 ** or less in temperature, and is ****** about all the processes so that viscous relaxation of an 
oxide film may not arise. And correlation with oxide film structure and Qbd which were predicted from 
the above-mentioned molecular dynamics calculation was investigated. 

[0452]Drawing 126 is the result of a Si-0 atom interval of the first contiguity searching for correlation 
with the frequency of occurrence of combination and Qbd used as 1.8 A or more, after computing a 
radial distribution function of a Si-0 atom of an oxide film in various stress fields. When the frequency 
of occurrence from which a Si-0 atom interval has become more than 1.8A exceeded 5% from Drawing 
126, it turned out that Qbd is decreasing notably and Si-zero-atom interval and an oxide film dielectric 

breakdown are related. The amounts of IR spectral shifts are about 15 cm -1 and intermediary ****. 
[0453]Change of amorphous NEFUTO work structure when compression stress is added to an oxide 
film, (1) Since density of an oxide film increases, two kinds of anticipation of anticipation that the degree 
of Si-O-Si bond angle decreases have been made from anticipation that Si-0 bond length decreases in 
number, IR measurement of (2) thin-film oxidizing film, and main KAMODERU. However, it newly found 
out that a rate of combination which increased density of an oxide film did not lead to Si-0 bond 
length's simple reduction, but bond length distribution spread rather, and was extended by a simulator 
which this invention persons developed uniquely as shown in Drawing 127 increased. In addition, it found 
out to correlation with Qbd. 

[0454]About a Reason for physical with correlation of Qbd and a postponed Si-0 combination, although 
this invention persons have not solved details yet, they guess - ** of the Reason as follows from a 
paper reported recently. A citation paper is "electron hole capture by local deformation and impurity in 
oxide film" p.101-104 reported by Mr. Kaneda with "a thin film and a surface physics subcommittee 
report of research, and formation, evaluation and reliability of ultra-thin silicon oxide." According to a 
place which Mr. Kaneda calculated using the molecular orbital method program GAUSSIAN "if there is 
extended Si-O-Si, capture of an electron hole to oxygen of the center will be induced. Since 



localization of the captured electron hole is strongly carried out into oxygen, to an electron, it works as 
a trapping center. As a result, many of oxygen captures an electron and they return to the original 
state. However, Si-0 combination which acquires several electron volts energy released by 
recombination of an electron and an electron hole, is extended and is weak is cut, and it jumps out to 
between lattices, for example, the part can generate a-like secondary defect like an oxygen deficiency. 
It is reported ". 

[0455]According to this invention persons, combination extended by stress addition increases like this 
report, and it is expected that an oxygen-deficiency defect is formed. Although various structures can 
be considered to an oxygen-deficiency defect, in structure where Si atom of oxygen-deficiency both 
ends was pulled apart greatly, an electronic state is also conjectured to change a lot, to form level in a 
band gap or to change to electron distribution which produces electrical conduction. And the insulation 
of an oxide film is imagined to deteriorate. 

[0456]Since an IR spectrum of an oxide film in inside of a stress field could be predicted and correlation 
with structure in an atom level and an electrical property was further acquired by making the above 
preparations, this invention persons tried "correction of a manufacturing process/' 
[0457]A Si wafer which progressed to a predetermined process was introduced into an oxide film 
forming device, and an oxide film of desired thickness was formed in 1050 ** and dry 0 2 atmosphere 

which contains chloride 5%. A wafer was quenched to a room temperature after oxide film formation, and 
an IR spectrum of a formed oxide film was measured. From measured value, a Si-0 asymmetric- 
stretching-vibration peak wave number of an IR spectrum was read, and when a shift amount from the 

wave number when stress has not arisen was calculated, a shift of 17 cm -1 was measured. 
[0458]From this shift amount, it was predicted that compression stress of about 0.2 GPa(s) is induced 

by oxide film in an oxide film. This shift amount was measured with shift critical-value [ of 15 cm ] ~ 1 . 
Like this heat process, since it was detected that it is over a critical value, a process was stopped. And 
stress relaxation by high temperature annealing was performed as correction of a manufacturing 
process. 

[0459]In order to compare an effect of this example, a high temperature annealing process was not 
added to some wafers, but a process was passed as it is. Annealing time, annealing temperature, a 
temperature falling speed, etc. which are the parameters of high temperature annealing for stress 
correction were computed in consideration of analysis of relaxation time by viscous relaxation prepared 
beforehand, and a re-diffusion action of an impurity diffused layer already formed into a substrate. In 
this case, 1000 ** and annealing for 30 seconds were performed. 

[0460]Since it was less than a critical value when an IR spectrum was measured again and a shift 
amount was computed after an end of annealing, a process was continued as it was. Qbd was measured 
by case where a process is corrected, and a case where it does not correct. As a result, in a case 
where it corrects, 45C/cm 2 and a case where it did not correct made 5 C/cm 2 , and an effect of 
process correction by this invention was acquired clearly. 

[0461 ]In a wafer which has another structure, gate oxide of desired thickness was formed in 950 ** dry 
0 2 atmosphere. When a wafer was quenched to a room temperature after oxide film formation and an IR 

spectrum of a formed oxide film was measured, a peak shift amount was 16 cm" 1 . Since this was over a 
critical value, an addition of annealing for stress relaxation was required for it, but when re-diffusion of 
an impurity diffused layer was taken into consideration, a solution of a suitable annealing condition did 
not exist. Then, it judged that process continuation was impossible for this wafer, and a wafer was 
discarded. 
[0462] 

[Effect of the Inventionjln the device, division single wafer processing, or the clustering manufacturing 
installation which it not only solves the problem of the monitoring facility in clustering, but manufactures 
semiconductor devices, such as MOSLSI, according to this invention as explained above, Even if the 
"spot" "spot" measured value of a limited quantity in each device is extended to the maximum and 
there is no test piece, it goes on, while a request corrects as a process. [ of measured value and limited 
qualitative contents ] According to this invention, fluctuation of (i) each process element, the check of 
whether to perform by diagnosing the superposition phenomenon of a gap of each element as the 



sequence of a (ii) request, and discovery of a (iii) change factor are also possible, furthermore — again 
— the feature of this system — change of a (iv) process — mathematical, the point it enabled it to deal 
with statistically and efficiently, and (v) — the point to which a new inference engine was made to add, 
and (vi) — it is in the point which added a new opposite direction system. 



[Translation done.] 
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[Brief Description of the Drawings] 

[Drawing 1]The figure showing the whole system by this invention. 

[Drawing 2]The figure showing that the inclination to each direction of a potential size serves as power 
of each direction. 

[Drawing 3]The figure showing the new abinitio molecular dynamics simulation system used for this 
invention. 

[Drawing 4]The figure showing the effect of starting an operation with the second atom (j) through the 
second G) from the third atom (k) further with first atomic (i). 

[Drawing 5]The figure making r the distance between particles and making theta a particle angle and in 
which showing the relation of the principal and accessory of these amounts of many. 
[Drawing 6]The figure showing positioning by the conventional calculation technique and the Takashina 
round partial differential method of this invention. 

[Drawing 7]The figure which explains cooperation with the process variation calculation unit I to be the 
molecular dynamics simulator II. 

[Drawing 8]The figure which explains cooperation with the process variation calculation unit I to be the 
molecular dynamics simulator II. 

[Drawing 9]The figure showing **** for an example of the input sequence of wide use 2 and a three- 
dimensional process simulator. 

[Drawing 10]The figure showing **** for an example of the input sequence of wide use 2 and a three- 
dimensional process simulator. 

[Drawing 1 1]The figure showing a part of binary tree by working example of this invention. 
[Drawing 12]The mimetic diagram of the Coulomb interaction at the time of ionization is shown. 
[Drawing 13]The mimetic diagram of the Coulomb interaction at the time of ionization is shown. 
[Drawing 14]The mimetic diagram of the Coulomb interaction at the time of ionization is shown. 
[Drawing 15]Easy atomic structure MODERU ** is shown. 
[Drawing 16]An alkaline metal halide molecule is shown. 
[Drawing 17]The table of the valence orbital index zeta is shown. 

[Drawing 18]The electric charge which dispelled and searched for the matrix (3063) is shown. 
[Drawing 19]The electric charge which dispelled and searched for the matrix (3063) is shown. 
[Drawing 20]The electric charge which dispelled and searched for the matrix (3063) is shown. 
[Drawing 21]The calculation result of Coulomb interaction J between Si-Si is shown. 
[Drawing 22]The electric charge of O is in 2 s, and the calculation result of JAB at the time of modeling 
the state where Si is in 2 s or 3 s is shown. 

[Drawing 23]What summarized the procedure of the newly developed above-mentioned electric charge 
numerical orientation method in the crystal system and the amorphous system is shown. 
[Drawing 24]The O-Si-0 angle which is an executed result of the molecular dynamics simulator by this 
invention, a Si-O-Si angle, and the data figure showing O-Si distance. 

[Drawing 25]The figure showing each subroutine including the variations calculation in the conventional 
simulation system. 



[Drawing 26]The figure showing the procedure of the variations calculation in the conventional 
simulation system. 

[Drawing 27]The characteristic figure showing the converging state of the conventional process and a 
change process. 

[Drawing 28]The data figure showing the accuracy of a variations. 
[Drawing 29]The data figure showing the accuracy of a variations. 
[Drawing 30]The data figure showing the accuracy of a variations. 

[Drawing 31]The figure showing the conventional impurity distribution by a simulation system and its 
example of change. 

[Drawing 32]The figure showing change and its execution combination of the process by the 
conventional simulation system. 

[Drawing 33]The figure showing some trees of search of the inference engine by the conventional 
simulation system. 

[Drawing 34]The figure of optimization of the perovskite by the new abinitioMD simulator which this 
invention persons developed, and Si system crystal. 

[Drawing 35]The construct figure of the simulator by original development by this invention person etc. 
[Drawing 36]The flow chart of the simulator by original development by this invention person etc. 
[Drawing 37]The figure of the example for C4v of the PTO perovskite by the simulator which this 
invention person etc. developed. 

[Drawing 38]The lineblock diagram of perovskite BaTi0 3 and PbTi0 3 . 

[Drawing 39]The potential output figure of the perovskite PTO by the simulator which this invention 
person etc. created newly. 

[Drawing 40]The example of creation of Si0 2 in the simulator which developed this invention person etc. 

[Drawing 41]The output figure of charge of the perovskite PTO by the simulator which this invention 
person etc. created newly. 

[Drawing 42]The example of distortion of the perovskite PTO by the abinitio/MD simulator which this 
invention person etc. created newly. 

[Drawing 43]The example of the oxygen deficiency of the perovskite PTO by the abinitio/MD simulator 
which this invention person etc. created newly. 

[Drawing 44]The example of movement of Pb of the perovskite PTO by the abinitio/MD simulator which 
this invention person etc. created newly. 

[Drawing 45]The example of movement of Pb of the perovskite PTO by the abinitio/MD simulator which 
this invention person etc. created newly. 

[Drawing 46]Some data of the simulator which this invention persons created 

[Drawing 47]The figure of saw YATAUWA. 

[Drawing 48]The figure of working example of this invention. 

[Drawing 49]The key map showing how the in-situ monitor measuring device is added in the clustering 
tool used for this invention. 

[Drawing 50]The figure showing signs that the method of making the limited measurand and sufficient 
information acquired from another side calculation contrasting is displayed as a picture [ graphics / on a 
monitor ]. 

[Drawing 51]The figure showing the temperature sequence as an example of the input data to the 
abinitio molecular dynamics simulation system by this invention. 

[Drawing 52]The key map showing execution prediction of the atom level by the abinitio molecular 
dynamics simulation system by this invention. 

[Drawing 53]The key map which visualized the atom level execution prediction by the abinitio molecular 
dynamics simulation system by this invention. 

[Drawing 54]The figure showing aging of the Si-O-Si angle which is an example of execution about Si0 2 
by the abinitio molecular dynamics simulation system by this invention. 

[Drawing 55]The prediction figure of measured value change of the temperature which is the example 
output to the input sequence of drawing 51. 

[Drawing 56]The distribution map of the coordination angle in membrane formation Si which is an 



example of atom level execution by the simulation system by this invention. 

[Drawing 57]The distribution map of the distance between Si-Si which is an example of atom level 
execution by the simulation system by this invention. 

[Drawing 58]The distribution map of the distance between Si-Si which is an example of atom level 
execution by the simulation system by this invention. 

[Drawing 59]The distribution map of the distance between Si-Si which is an example of atom level 
execution by the simulation system by this invention. 

[Drawing 60]The distribution map of Si coordination angle which is an example of atom level execution 
by the simulation system by this invention. 

[Drawing 61]The distribution map of the distance between Si-Si which is an example of atom level 
execution by the simulation system by this invention. 

[Drawing 62]The example of atom level execution by the simulation system by this invention. The 
distribution map of the distance between Si-Si. 

[Drawing 63]The distribution map of Si coordination angle which is an example of atom level execution 
by the simulation system by this invention. 

[Drawing 64]The distribution map of Si coordination angle which is an example of atom level execution 
by the simulation system by this invention. 

[Drawing 65]The distribution map of Si coordination angle which is an example of atom level execution 
by the simulation system by this invention. 

[Drawing 66]The distribution map of the distance between Si-Si which is an example of atom level 
execution by the simulation system by this invention. 

[Drawing 67]The distribution map of the distance between Si-Si which is an example of atom level 
execution by the simulation system by this invention. 

[Drawing 68]The distribution map of the distance between Si-Si which is an example of atom level 
execution by the simulation system by this invention. 

[Drawing 69]The distribution map of the distance between Si-Si which is an example of atom level 
execution by the simulation system by this invention. 

[Drawing 70]The figure showing the another whole system constituted by this invention. 
[Drawing 71]The characteristic figure showing the conventional Raman measurement result. 
[Drawing 72]The characteristic figure showing the conventional Raman measurement result. 
[Drawing 73]The characteristic figure showing the conventional Raman measurement result. 
[Drawing 74]The figure showing the model for explaining the conventional Raman measurement result. 
[Drawing 75]The characteristic figure showing the strength reduction produced by the characteristic 
figure and stress which show the wavelength of normal mode vibration of the Si0 4 tetrahedron as a 

function of an O-Si-0 bond angle. 

[Drawing 76]The model figure showing the defect in Si0 2 . 

[Drawing 77]The characteristic figure showing the example of execution of the molecular dynamics 
simulator by this invention. 

[Drawing 78]The characteristic figure showing the example of execution of the molecular dynamics 
simulator by this invention. 

[Drawing 79]The characteristic figure showing the example of execution of the molecular dynamics 
simulator by this invention. 

[Drawing 80]The figure showing a motion of the atom which is an example of execution of the molecular 
dynamics simulator by this invention. 

[Drawing 81]The characteristic figure showing the X-ray-spectrum prediction result which is an 
example of execution of the molecular dynamics simulator by this invention. 

[Drawing 82]The characteristic figure showing the stress prediction result which is an example of 
execution of the molecular dynamics simulator by this invention. 

[Drawing 83]The characteristic figure showing the computed result by the molecular dynamics simulator 
which this invention persons developed newly. 

[Drawing 84]The figure showing the consistent calculation including change of the molecular dynamics 
simulator using the variations technique by this invention. 

[Drawing 85]The figure showing the search technique by the molecular dynamics simulator by this 



invention. 

[Drawing 86]The characteristic figure showing the situation of the temperature setting at the time of 
the amorphous Si0 2 creation by this invention persons. 

[Drawing 87]The characteristic figure showing the pressure variation at the time of amorphous Si0 2 
creation in the computer by this invention persons. 

[Drawing 88]The characteristic figure showing the volume change at the time of the amorphous Si0 2 

creation by the molecular dynamics simulator which this invention persons created. 

[Drawing 89]The characteristic figure showing the volume change at the time of the amorphous Si0 2 

creation by the molecular dynamics simulator which this invention persons developed newly. 

[Drawing 90]The characteristic figure showing the potential energy change at the time of the amorphous 

Si0 2 creation by the molecular dynamics simulator which this invention persons developed newly. 

[Drawing 91]The characteristic figure showing the temporal change of two sides in the calculation cells 
at the time of the amorphous Si0 2 creation by the molecular dynamics simulator which this invention 

persons developed newly. 

[Drawing 92]The characteristic figure showing change of the calculated energy for the molecular 
dynamics simulator which this invention persons developed newly using the algorithm of the usual 
Verlet. 

[Drawing 93]The characteristic figure using the algorithm of advanced Verlet for the molecular dynamics 
simulator which this invention persons developed newly, and showing change of the calculated energy 
for it. 

[Drawing 94]The characteristic figure showing the temporal change of the dipole moment calculated 
using the molecular dynamics simulator which this invention persons developed newly. 
[Drawing 95]The characteristic figure showing the power spectrum for which this invention persons 
asked further from the dipole moment using the molecular dynamics simulator developed newly. 
[Drawing 96]The figure applying the molecular dynamics simulator which this invention persons 
developed newly to alpha-SiO 2 and in which showing the calculated IR spectrum. 

[Drawing 97]The figure showing the IR spectrum as which it was predicted when the molecular dynamics 
simulator which this invention persons developed newly was applied to alpha-SiO 2 and the external 
pressure was added. 

[Drawing 98]The characteristic figure showing the volume change at the time of amorphous SiO by the 
molecular dynamics simulator which this invention persons developed newly, and creation. 
[Drawing 99]The characteristic figure showing change of the calculated energy for the molecular 
dynamics simulator which this invention persons developed newly using the algorithm of the usual 
Verlet. 

[Drawing 100]The figure showing the computational procedure by this invention persons. 
[Drawing 101]The figure showing the computational procedure by this invention persons. 
[Drawing 102]The figure showing the sliding surface in Si (100), decomposition shearing stress, etc. 
[Drawing 103]The mimetic diagram showing a reaction near the Si0 2 and Si/Si0 2 interface. 

[Drawing 104]The figure showing serial calculation of viscoelasticity calculation typically. 

[Drawing 105]The temperature of the thickness of an oxide film, the characteristic figure showing time 

dependency. 

[Drawing 106]The temperature of the thickness of an oxide film, the characteristic figure showing time 
dependency. 

[Drawing 107]The temperature of the thickness of an oxide film, the characteristic figure showing time 
dependency. 

[Drawing 108]The figure showing the principal stress distribution and maximum shear stress in trench 
oxidation at the time of using the doubled coefficient. 

[Drawing 109]The characteristic figure showing change of the temperature which can be set like a 
typical heat process, and accumulation stress. 

[Drawing 110]The calculation prediction figure showing distribution of sigmaxx in a trench angle. 
[Drawing 1 1 1]The calculation prediction figure showing distribution of sigmaxx in a trench angle. 



[Drawing 112]The calculation prediction figure showing distribution of sigmaxx in a trench angle. 
[Drawing 1 13]The calculation prediction figure showing distribution of sigmaxx in a trench angle. 
[Drawing 114]The calculation prediction figure showing distribution of sigmaxx in a trench angle. 
[Drawing 1 15]The calculation prediction figure showing distribution of sigmaxx in a trench angle. 
[Drawing 1 16]The figure showing the shape of the optimized oxidation prediction. 
[Drawing 1 1 7]The sectional view showing the formation process of the gate oxidizing film in working 
example 7. 

[Drawing 118]The sectional view showing the formation process of the gate oxidizing film in working 
example 7. 

[Drawing 1 19]The characteristic figure showing the temperature rising step after the gate oxidizing film 
formation in working example 7. 

[Drawing 120]The characteristic figure showing change of the stress by difference of a temperature 
rising step 

[Drawing 121]The characteristic figure showing change of the stress by difference of a temperature 
rising step. 

[Drawing 122]The figure showing the example computation of the stress distribution. 

[Drawing 123]The figure showing the relation between the oxidizing temperature in working example 7, 

and the compression stress in a film. 

[Drawing 124]The figure showing the calculation result of change of the IR spectrum by stress addition. 
[Drawing 125]The figure showing the relation between the compression stress in an oxide film, and the 
wave number shift of a Si-O-Si asymmetric stretch peak. 

[Drawing 126]The figure in which an about [ 1st ] Si-0 atom interval shows the relation of the 
frequency of occurrence of combination and Qbd used as 1.8 A or more. 

[Drawing 127]The figure showing change of the Si-0 atom interval distribution by the existence of 
stress addition. 

[Drawing 128]The figure showing the relation between the generation of dRAM, and chip cost. 
[Drawing 129]The mimetic diagram showing the manufacturing process of the semiconductor device by 
a conventional example. 

[Drawing 130]The figure showing the clustering tool by a conventional example. 

[Drawing 131]The figure showing the composition of the information processing system which realizes 
the process simulation by this invention. 

[Drawing 132]The figure showing the example of a three-dimensional rigid matrix. 



[Translation done.] 
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